!     Last change:  SYR   3 Sep 2003    2:41 pm
!==============================================================================
!          PROGRAMME: COLBERT
!==============================================================================

!==============================================================================

!     ------------------------------------------------------------------
!     ------------------------------------------------------------------
!     T A B L E   D E S   M A T I E R E S
!     ------------------------------------------------------------------
!     ------------------------------------------------------------------

!     INITIALISATIONS
!     DEBUT
!     SETTING
!     COMP(ifi,iarg)
!     DECOMP(iarg,ifi)
!     TASSE(ilit,itas)
!     DETASSE(itas,ilit)

!     PROGRAMMES DE CALCUL LITTERAL
!     RAZX
!     RAZ(n)
!     TRIX
!     TRIR
!     TRIARG
!     DECAL(long)
!     TRIFI(n1,n2)
!     DECFI(n,long)
!     PURGE(n)
!     MULSER(n1,n2,n3,natt,factr)
!     SOMSER(n1,n2,natt,fact)
!     SOMCTE(n1,n2,natt,fact)
!     TRANS(n1,n2)
!     DERIVE(n1,n2)
!     PRIMIT(n1,n2)
!     DERARG(n1,n2)
!     DERLIT(n1,n2)
!     SUBTIME(n1,n2)

!     IMPRESSION
!     IMPRFLOT(n)
!     IMPRFRAC(n)
!     FRACONT(f,inum,iden,fail)

!     SERIES DU PROBLEME DES 2-CORPS
!     INIHANS
!     HANSEN(n,m,k,ordre,h)  [E]
!     DEVEL(na,imax,f)
!     RAC(a,imax,f)
!     HYPERG(a,b,c,imax,f)
!     DEVBET(imax,f)
!     FACT(k,n)  [F]
!     COEF_P(n,s,m,nu,imax,p)
!     PROD(p1,p2,imax,p)
!     DEUXCORPS(n,m,xcos,ordre,iser)

!     ENTREE ET SORTIE DE SERIES
!     SINGLE(n,m,xcos,iser)
!     SAISIE_F(nomdata,iser)
!     SORTIE_F(nomdata,n)
!     OUTGREG(n1,n2)

!     INTERPRETEUR DE COMMANDES (DECODAGE DES INSTRUCTIONS)
!     INTERPRETE
!     ANAGAUCHE(instruc,xfor)
!     ANADROITE(instruc)
!     CYCLOPE(ncycle,ntype,xout)
!     SYMBOLE(chaine,nom,nx)
!     NOMBRE(chaine,dx,nx)
!     EXPLORE(instruc)
!     SICONST(chaine)

!     GESTION DES TABLEAUX DE NOMS ET DE NUMEROS
!     INITAB
!     TRIN(nom)
!     KILLN(nom)
!     AJOUTN(nom)
!     PUSH(n)
!     POP(n)
!     TRIC(nom)
!     KILLC(nom)
!     AJOUTC(nom,ct)
!     QUELNOM(n,nom)

!     GESTION DISQUE
!     CASH(nin,xlec,nout)
!     EDISK(n,ndisk)
!     LDISK(n,ndisk)
!
!     GESTION DES ERREURS
!     ERREUR(n)

!     Liste des routines utlisateur (30 - 07 - 03)
!        CHERCHEC
!        CONSER
!        HANSENC
!        HANSENS
!        SINGLEC
!        SINGLES
!        DERARG
!        DERLIT
!==============================================================================

      MODULE PARTLIT

!     Parametrage des listes sur la ligne de commande
!     ICOMMANDE=1 : FORTRAN LAHEY ; ICOMMANDE=2 : FORTRAN HP UNIX
      integer, parameter      :: ICOMMANDE=1                            !$$$$$

!----------------------------------------------------------------------
!     Unites logiques
      integer, parameter      :: NLOG1=3   ! Moniteur (sequentiel) : 
!                                            'MONITEST'
      integer, parameter      :: NLOG2=4   ! Fichier de saisie 
!                                            en lecture|ecriture (sequentiel)
      integer, parameter      :: NLOG3=15  ! Fichier cash (direct) : 
!                                            'TAMPON'
      integer, parameter      :: NLOG6=10  ! Fichier des parametres utiles : 
!                                            'SETTING.DAT'
      integer, parameter      :: NLOG10=30 ! Routines  (sequentiel) : 
!                                            'MONITEST.SUB'
      integer                 :: NLOGY     ! Copie unite du moniteur|routines

      integer, parameter      :: MXSER=128
      integer, parameter      :: NBT=32
      integer, parameter      :: MXTAMP=MXSER-NBT
      integer, parameter      :: MXCTE=128

      integer, parameter      :: MXDIMR=4096
      integer, parameter      :: MXDIMT=4096
      integer, parameter      :: MXDIMX=32

      integer, parameter      :: MXNLIT=10
      integer, parameter      :: MXNARG=10
      integer                 :: NLIT,NARG

      integer, parameter      :: ISTOP=16
      integer, parameter      :: NPUIS2(ISTOP)=(/
     & 16384,8192,4096,2048,1024,512,256,128,64,32,16,8,4,2,1,0/)
      character*20, parameter :: BIGNOM=''

      integer, parameter      :: MAXINT=2147483647

!     Tableaux R, T et X
      integer, dimension(MXDIMR)          :: lr
      integer, dimension(MXDIMR)          :: jr
      real*8, dimension(MXDIMR)           :: r

      integer, dimension(MXDIMT,NBT)      :: lt
      integer, dimension(MXDIMT,NBT)      :: jt
      real*8, dimension(MXDIMT,NBT)       :: t
      integer, dimension(NBT)             :: ntt

      integer, dimension(MXDIMX)          :: lx
      integer, dimension(MXDIMX)          :: jx
      real*8, dimension(MXDIMX)           :: x
      integer                             :: ntx

      logical*1                     :: trigr
      integer                       :: kdimr,milr,idebr,ntr,nr,irangr
      logical*1, dimension(NBT)     :: trig
      integer                       :: kdimt,milt,idebt

      integer, dimension(MXSER)           :: numser
      character*20, dimension(MXSER)      :: nomser
      integer                             :: kdimser,milser,idebser
      integer                             :: ntser,nser
      logical*1                           :: xser

      real*8, dimension(MXCTE)            :: cte
      character*20, dimension(MXCTE)      :: nomcte
      integer                             :: kdimcte,milcte,idebcte
      integer                             :: ntcte,ncte
      logical*1                           :: xcte,xtot

      integer, dimension(NBT)             :: ndismem,nmem
      integer, dimension(MXTAMP)          :: ndistamp
      integer                             :: ntopmem,ntoptamp,ntmem

      integer                             :: lw,jw
      real*8                              :: rw

      logical*1                           :: sifrac,siargu,sishort
      integer                             :: maxout,maximpr
      real*8                              :: prec

      integer, dimension(0:MXNLIT)        :: mulit,mxpar,ipds
      integer, dimension(MXNARG)          :: mularg
      integer                             :: miarg,mzer
      integer                             :: idiv,maxord,maxop

      integer                             :: ipl,ipf

      character*101                       :: zrecbis
      character*100                       :: instrucbis

      character*20                        :: gauche
      logical*1                           :: xcteg,xserg
      integer                             :: ncop
      integer, dimension(0:2)             :: nop
      real*8                              :: cop

      logical*1                           :: xfor, xiter
      integer                             :: ITER, LIGNE

      character*12                        :: MONIT, FILSUB
      character*20                        :: subutil

      character*20                        :: titre
      character*50                        :: zimp,zimp1,zimp2,zimp3
!-----------------------------------------------------------------------
!     Les phases PHASE() et les frequences FREQ()
!     sont chargees dans la routine SETTING a partir des donnees du fichier
!     'SETTING.DAT'
!     Ces donnees definissent les unites. Les angles sont en radian
!     les frequence sont en radian/unite de temps (annee, siecle, millenaire...)
!---------------------------------------------------------------------

      real*8, dimension(MXNARG)           :: FREQ, PHASE
      real*8, dimension(MXNLIT)           :: VAR

      END MODULE PARTLIT

!==============================================================================

      PROGRAM COLBERT
      use PARTLIT
      implicit real*8 (a-h,o-z)
!**********************************************************************

!  ------------------ M O N I T E U R ------------------

!**********************************************************************

!----------------------------------------------------------------------
!     Manipulateur de series de Fourier a coefficients litteraux
!----------------------------------------------------------------------

      call SETTING
      call INITAB
      call DEBUT
      call INIHANS

!     ******************************************************************
      prec=1.d-10
      maximpr=4
      sifrac=.TRUE.
      maxop=maxord

      call INTERPRETE
!     print *,'Bonjour...'
!     call INTERPRETE
!     print *,'Bonsoir...'
      stop 'Termine...'

      END PROGRAM COLBERT



!***********************************************************************
!***********************************************************************

!     --------------- I N I T I A L I S A T I O N S ----------------

!***********************************************************************
!***********************************************************************

!----------------------------------------------------------------------

!                          D E B U T

!----------------------------------------------------------------------
!     Demarrage des sous-programmes
!----------------------------------------------------------------------

      subroutine DEBUT
      use PARTLIT
      implicit real*8 (a-h,o-z)
!----------------------------------------------------------------------

      character*3          :: cw
      character*12         :: cfile
      character*12         :: TAMPON
      character*80         :: liste
      character*20         :: word

      character*20         :: z1="(1xi4,1x,'[',i2,']',"
      character*4          :: z13="(1X,"

      character*17         :: z2="i3,(1x,d15.8),2x"
      character*17         :: z21="i3,(3x,d17.10),2x"
      character*17         :: z22="i3,i10,'/',i9,2x"
      character*16         :: z23="i3,1x,d23.16,2x,"

      character*3          :: z3="i3)"
      character*12         :: z4="x,a3)"

      integer, parameter   :: its=48
!----------------------------------------------------------------------
!     Definition de INFR et INFT, rangs dans le tableau PUIS2(),
!     qui donnent les dimensions actuelles de R() <= MXDIMR et
!     T() <= MXDIMT respectivement.
!----------------------------------------------------------------------
      INFTAB(k)=15-LOG(REAL(k,KIND=8))/LOG(2.D0)
      infr=inftab(MXDIMR)
      inft=inftab(MXDIMT)
      kdimr=NPUIS2(infr)
      kdimt=NPUIS2(inft)
      milr=kdimr/2
      milt=kdimt/2
      idebr=infr+2
      idebt=inft+2
!----------------------------------------------------------------------
!     Initialisation des parametres utiles a la manipulation des
!     tables des noms et ces constantes.
!----------------------------------------------------------------------
      infz=inftab(MXSER)
      kdimser=NPUIS2(infz)
      milser=kdimser/2
      idebser=infz+2
      infz=inftab(MXCTE)
      kdimcte=NPUIS2(infz)
      milcte=kdimcte/2
      idebcte=infz+2

      ipl=1          ! valeur par defaut
      ipf=1             ! valeur par defaut

      if (NARG > 10) call ERREUR(200)
      if (NARG > MXNARG) call ERREUR(200)
      if (NLIT > 10) call ERREUR(201)
      if (NLIT > MXNLIT) call ERREUR(201)

! Initialisation des arguments
!----------------------------------------------------------------------
!     Initialisation pour la compression-decompression d'arguments
!     sur un entier (INTEGER*4). Le nombre de composantes de l'argument
!     est fixe par NARG.
!----------------------------------------------------------------------
      if (NARG == 1) then
         miarg=MAXINT
      else
         a=LOG(DFLOAT(MAXINT))
         iw=EXP(a/NARG)
         miarg=(iw-1)/2
         mbaz=2*miarg+1
         mularg(NARG)=1
         mzer=1
         do i=NARG-1,1,-1
         mularg(i)=mularg(i+1)*mbaz
         mzer=mzer+mularg(i)
         end do
         mzer=miarg*mzer
      end if
      write (*,100) NARG,miarg
100   format(' Nombre de composantes',
     + ' de l''argument (NARG) : ',i2,/
     + ' --> Chacune des composantes de l''argument doit',
     + ' etre <=',i10,' en valeur absolue.')

!----------------------------------------------------------------------
!     Definition des parametres lies a la valeur de MAXORD :
!     Initialisation des parametres litteraux (n-uplets) admissibles
!     sur un entier. Leur nombre est donne par NLIT. L'ordre maximum
!     est fixe par MAXORD; les poids associes aux parametres sont lus
!     dans IPDS(*). La routine calcule IDIV; IDIV est le diviseur entier
!     d'un ensemble de parametres litteraux (ou n-uplets) compactes
!     sur un entier (i.e. ITAS) permettant le calcul de son ordre
!     avec: ITAS/NDIV.
!----------------------------------------------------------------------
      do i=1,NLIT
      mxpar(i)=maxord/ipds(i)
      end do
      write (*,200) NLIT,maxord
200   format(' Nombre de variables litterales',
     + ' (NLIT) : ' ,i2,' ; Ordre maximum (MAXORD) : ',i2)
      write (*,201) (mxpar(i), i=1,NLIT)
201   format (' Parametres:',10i4)

      supint=MAXINT
      w=1
      mulit(NLIT)=1
      do i=NLIT,1,-1
      w=w*(mxpar(i)+1)
      mulit(i-1)=w
      if (w > supint) call ERREUR(202)
      end do
      idiv=mulit(0)

!----------------------------------------------------------------------
      zimp=z1//char(its+NLIT)//z2//char(its+NARG)//z3
      zimp1=z1//char(its+NLIT)//z21//char(its+NARG)//z3
      zimp2=z1//char(its+NLIT)//z22//char(its+NARG)//z3
      zimp3=z13//char(its+NLIT)//z23//char(its+NARG)//z3
      nb=14+3*NLIT
      n1=nb/10
      n2=nb-10*n1
      titre="("//char(its+n1)//char(its+n2)//z4
!----------------------------------------------------------------------

      NLOGY=NLOG1

      MONIT='monitest'                                                  !$$$$$
      FILSUB='monitest.sub'                                             !$$$$$
      TAMPON='tampon'                                                   !$$$$$
      sishort=.FALSE.



!-----------------------------------------------------------------------
!     ATTENTION:
!     La routine GETARG reagit de maniere  differente selon
!     le compilateur.
!     Avec le compilateur LAHEY, la liste des arguments qui suivent
!     la commande commence au numero 1.
!     Sous UNIX la liste commence au numero 2
!-----------------------------------------------------------------------
      liste=''
      do i=ICOMMANDE,ICOMMANDE+4   ! Liste des parametres de GETARG

         call GETARG(i,word)

         n=LEN_TRIM(liste)

          if (i==ICOMMANDE) then
            liste=word
         else
            liste=liste(1:n)//' '//word
         end if

      end do

!---------------------------------------------------------------------

      m=INDEX(liste,'-t')
      if (m /=0) then
         STOP 'L''option trace -t n''existe pas dans COLBERT'
      end if
      liste=ADJUSTL(liste)
      n=INDEX(liste,' ')
      if (n == 1) goto 10
      if (liste(1:1) /= '-') then
         MONIT=liste(1:n-1)
         liste=ADJUSTL(liste(n:))
         n=INDEX(liste,' ')
         if (n == 1) goto 10
      end if
      do i=1,2
         cw=liste(1:3)
         liste=liste(4:)
         liste=ADJUSTL(liste)
         n=INDEX(liste,' ')
         if (n == 1) exit
         cfile=liste(:n)
         if (cw == '-r:')
     & STOP 'La sauvegarge (IN, option -r) n''existe pas dans COLBERT'
         if (cw == '-w:')
     & STOP 'La sauvegarge (OUT, option -w) n'' existe pas dans COLBERT'
         liste=ADJUSTL(liste(n:))
      end do
10    continue
      FILSUB=TRIM(MONIT)//'.SUB'

      OPEN (NLOG1,FILE=MONIT)

      IBYTE=8*MXDIMT
      OPEN (NLOG3,access='DIRECT',file='TAMPON',recl=IBYTE)

      END subroutine DEBUT

!---------------------------------------------------------------------

!                             S E T T I N G

!---------------------------------------------------------------------
!     Chargement des tableaux des frequences FREQ(), des phases PHASE()
!-----------------------------------------------------------------------
      subroutine SETTING

      use PARTLIT
      implicit real*8      (a-h,o-z)
!---------------------------------------------------------------------

      NAMELIST /FP/ NLIT,NARG,MAXORD,IPDS,FREQ,PHASE,VAR

      OPEN (UNIT=NLOG6,FILE='SETTING.DAT')
      read (NLOG6,NML=FP)

      END subroutine SETTING

!----------------------------------------------------------------------

!                          C O M P

!----------------------------------------------------------------------
!     Compression d'un argument a NARG<=10 composantes, sur un entier.
!     Entree: IFI(NARG), argument en clair; Sortie: IARG, argument
!     comprime sous forme interne (entier).
!----------------------------------------------------------------------
      subroutine COMP(ifi,iarg)
      use PARTLIT
      implicit real*8 (a-h,o-z)

      integer, dimension(NARG)            :: ifi
!----------------------------------------------------------------------
      if (NARG == 1) then
      iarg=ifi(1)
      RETURN
      end if
      iarg=0
      do i=1,NARG
      iarg=iarg+ifi(i)*mularg(i)
      end do
      RETURN
      END subroutine COMP

!----------------------------------------------------------------------

!                          D E C O M P

!----------------------------------------------------------------------
!     Decompression d'un argument a NARG composantes, comprime sur un
!     entier. Entree: IARG, argument comprime sous forme interne;
!     Sortie: IFI(NARG), argument en clair decomprime.
!----------------------------------------------------------------------
      subroutine DECOMP(iarg,ifi)
      use PARTLIT
      implicit real*8 (a-h,o-z)

      integer, dimension(NARG)            :: ifi
!----------------------------------------------------------------------

      if (NARG == 1) then
      ifi(1)=iarg
      RETURN
      end if
      jarg=iarg+mzer
      do i=1,NARG
      ifi(i)=jarg/mularg(i)-miarg
      jarg=jarg-(ifi(i)+miarg)*mularg(i)
      end do
      RETURN
      END subroutine DECOMP

!----------------------------------------------------------------------

!                          T A S S E

!----------------------------------------------------------------------
!     Compression d'un argument litteral a NLIT composantes, sur un
!     entier. Entree: ILIT(1:NLIT), argument en clair;
!     Sortie: ILIT(0), ordre calcule avec les poids associes;
!     ITAS, argument comprime sous forme interne.
!----------------------------------------------------------------------
      subroutine TASSE(ilit,itas)
      use PARTLIT
      implicit real*8 (a-h,o-z)

      integer, dimension(0:NLIT)          :: ilit

!----------------------------------------------------------------------
      ilit(0)=0
      do i=1,NLIT
      ilit(0)=ilit(0)+ipds(i)*ilit(i)
      end do
      if (ilit(0) > maxop) call ERREUR(203)
      itas=ilit(NLIT)
      do i=0,NLIT-1
      itas=itas+ilit(i)*mulit(i)
      end do
      RETURN
      END subroutine TASSE

!----------------------------------------------------------------------

!                          D E T A S S E

!----------------------------------------------------------------------
!     Decompression d'un argument a NLIT composantes, comprime sur un
!     entier. Entree: ITAS, argument comprime sous forme interne;
!     Sortie: ILIT(0:NLIT), argument litteral en clair decomprime.
!----------------------------------------------------------------------
      subroutine DETASSE(itas,ilit)
      use PARTLIT
      implicit real*8 (a-h,o-z)

      integer, dimension(0:NLIT)          :: ilit

!----------------------------------------------------------------------
      jtas=itas
      do i=0,NLIT
      ilit(i)=jtas/mulit(i)
      jtas=jtas-ilit(i)*mulit(i)
      end do
      if (ilit(0) > maxop) call ERREUR(204)
      RETURN
      END subroutine DETASSE

!***********************************************************************
!***********************************************************************

!     -- P R O G R A M M E S  D E  C A L C U L  L I T T E R A L --

!***********************************************************************
!***********************************************************************

!----------------------------------------------------------------------

!                          R A Z X

!----------------------------------------------------------------------
!     Remise a zero de X(*), de la table des arguments JX(*) a MAXINT
!     et des n-uplets LX(*) a MAXINT.
!----------------------------------------------------------------------
      subroutine RAZX
      use PARTLIT
      implicit real*8 (a-h,o-z)

!----------------------------------------------------------------------
      ntx=0
      do 1 i=1,32
      x(i)=0
      lx(i)=MAXINT
1     jx(i)=MAXINT
      RETURN
      END subroutine RAZX

!----------------------------------------------------------------------

!                          R A Z

!----------------------------------------------------------------------
!     Remise a zero d'un tableau T(*,N) d'indice N si N <>0.
!     Remise a zero de R(*) si N=0.
!----------------------------------------------------------------------
      subroutine RAZ(n)
      use PARTLIT
      implicit real*8 (a-h,o-z)

!----------------------------------------------------------------------
      if (n == 0) then
      ntr=0
      do 1 i=1,kdimr
         r(i)=0
         lr(i)=MAXINT
1        jr(i)=MAXINT
      RETURN
      endif
      ntt(n)=0
      do 2 i=1,kdimt
       t(i,n)=0
       lt(i,n)=MAXINT
2      jt(i,n)=MAXINT
      RETURN
      END subroutine RAZ

!----------------------------------------------------------------------

!                          T R I X

!----------------------------------------------------------------------
!     Tri dans la table X(*) par n-uplets puis par arguments croissants.
!     Entree: terme (RW), argument (JW), n-uplet (LW) a ranger dans X(*).
!     Sortie: si JW et LW existent dans X(*) le terme est ajoute,
!     sinon il est intercale. Le rangement est opere d'abord sur les
!     n-uplets croissants, ensuite sur les arguments croissants.
!----------------------------------------------------------------------
      subroutine TRIX
      use PARTLIT
      implicit real*8 (a-h,o-z)

!----------------------------------------------------------------------
      nx=16
      do 1 i=12,16
      if (jw == jx(nx).and.lw == lx(nx)) then
         x(nx)=x(nx)+rw
         RETURN
      endif
      isix=lx(nx)
      if (lw > isix) then
         nx=nx+NPUIS2(i)
         l=1
      else
         if (lw == isix) then
            if (jw > jx(nx)) then
            nx=nx+NPUIS2(i)
            l=1
            else
            nx=nx-NPUIS2(i)
            l=0
            endif
         else
            nx=nx-NPUIS2(i)
            l=0
         endif
      endif
1     continue
      nx=nx+l
      ntx=ntx+1
      if (nx == ntx) then
      lx(nx)=lw
      jx(nx)=jw
      x(nx)=rw
      RETURN
      endif
      imax=ntx-nx
      if (imax <= 0) call ERREUR(205)
      do 2 i=1,imax
      ib=ntx-i
      ia=ib+1
      lx(ia)=lx(ib)
      jx(ia)=jx(ib)
 2    x(ia)=x(ib)
      lx(nx)=lw
      jx(nx)=jw
      x(nx)=rw
      RETURN
      END subroutine TRIX

!----------------------------------------------------------------------

!                          T R I R

!----------------------------------------------------------------------
!     Tri dans la table R(*) par n-uplets puis par arguments croissants.
!     Entree: terme (RW), argument (JW), n-uplet (LW) a ranger dans R(*);
!     Sortie: si JW et LW existent dans R(*) (ou JW > JR max.) alors
!     IRANGR=1, le terme est ajoute (NTR=NTR+1). Sinon, JW et/ou LW
!     n'existent pas, et IRANGR=0; il n'y a pas d'interclassement.
!     NR est le rang du terme dans R(*); NTR est inchange.
!----------------------------------------------------------------------
      subroutine TRIR
      use PARTLIT
      implicit real*8 (a-h,o-z)

!----------------------------------------------------------------------
      nr=milr
      irangr=1
      do 1 i=idebr,16
      if (jw == jr(nr).and.lw == lr(nr)) then
         r(nr)=r(nr)+rw
         RETURN
      endif
      isr=lr(nr)
      if (lw > isr) then
         nr=nr+NPUIS2(i)
         l=1
      else
         if(lw == isr) then
            if (jw > jr(nr)) then
            nr=nr+NPUIS2(i)
            l=1
            else
            nr=nr-NPUIS2(i)
            l=0
            endif
         else
            nr=nr-NPUIS2(i)
            l=0
         endif
      endif
1     continue
      nr=nr+l
      ntr=ntr+1
      if (nr == ntr) then
      lr(nr)=lw
      jr(nr)=jw
      r(nr)=rw
      RETURN
      endif
      ntr=ntr-1
      irangr=0
      RETURN
      END subroutine TRIR

!----------------------------------------------------------------------

!                          T R I A R G

!----------------------------------------------------------------------
!     Tri d'un n-uplet (LW) et d'un argument (JW) quelconque;
!     (RW) est le coefficient associe. Le terme est ajoute a R(*)
!     s'il existe, sinon il est ajoute a X(*) et/ou interclasse.
!----------------------------------------------------------------------
      subroutine TRIARG
      use PARTLIT
      implicit real*8 (a-h,o-z)

!----------------------------------------------------------------------
      call TRIR
      if (irangr == 0) call TRIX
      if ((ntr+ntx) >= kdimr) call ERREUR(206)
      if (ntx == 32) call DECAL (32)
      RETURN
      END subroutine TRIARG

!----------------------------------------------------------------------

!                          D E C A L

!----------------------------------------------------------------------
!     Interclassement des elements de X(*) dans R(*); LONG: nombre
!     d'elements a classer.
!----------------------------------------------------------------------
      subroutine DECAL(long)
      use PARTLIT
      implicit real*8 (a-h,o-z)

!----------------------------------------------------------------------
      if (long == 0) RETURN
      m=ntr
      do 1 i=0,long-1
      k=long-i
      rw=x(k)
      jw=jx(k)
      lw=lx(k)
      call TRIR
      m=m+1
      jmax=m-nr
      if (jmax <= 0) call ERREUR(207)
      do 2 j=1,jmax
         ib=m-j
         ia=ib+k
         lr(ia)=lr(ib)
         jr(ia)=jr(ib)
2        r(ia)=r(ib)
      lr(nr)=lw
      jr(nr)=jw
      r(nr)=rw
1     m=nr
      ntr=ntr+long
      call RAZX
      RETURN
      END subroutine DECAL

!----------------------------------------------------------------------

!                          T R I F I

!----------------------------------------------------------------------
!     Tri sur les arguments et les n-uplets dans le tableau N1 (ou 0);
!     le resultat est range dans le tableau de numero N2.
!     Le tri est effectue d'abord sur l'ordre croissant des arguments,
!     ensuite sur l'ordre croissant des n-uplets.
!     X(*) sert de tableau de manoeuvre intermediaire.
!     ATTENTION: Ce tri n'a d'utilite que pour ordonner les series pour
!     l'impression par arguments croissants. De telles serie triees
!     ne peuvent pas etre utilisees pour des operations arithmetiques
!     telles que MULSER ou SOMSER.
!----------------------------------------------------------------------
      subroutine TRIFI(nn1,nn2)
      use PARTLIT
      implicit real*8 (a-h,o-z)

!----------------------------------------------------------------------
      call CASH(nn1,.TRUE.,n1,0)
      call CASH(nn2,.FALSE.,n2,0)
      call RAZX
      call RAZ(n2)
      if (n1 == 0) then
      jmax=ntr
      trig(n2)=trigr
      else
      jmax=ntt(n1)
      trig(n2)=trig(n1)
      endif
      do j=1,jmax
      if (n1 == 0) then
         jw=jr(j)
         lw=lr(j)
         rw=r(j)
      else
         jw=jt(j,n1)
         lw=lt(j,n1)
         rw=t(j,n1)
      endif
      nx=16
      do i=12,16
         isix=jx(nx)
         if (jw > isix) then
            nx=nx+NPUIS2(i)
            l=1
         else
            if (jw == isix) then
            if (lw > lx(nx)) then
               nx=nx+NPUIS2(i)
               l=1
            else
               nx=nx-NPUIS2(i)
               l=0
            endif
            else
            nx=nx-NPUIS2(i)
            l=0
            endif
         endif
      end do
      nx=nx+l
      ntx=ntx+1
      if (nx <  ntx) then
         imax=ntx-nx
         do i=1,imax
            ib=ntx-i
            ia=ib+1
            lx(ia)=lx(ib)
            jx(ia)=jx(ib)
            x(ia)=x(ib)
         end do
      else
         if (nx > ntx) call ERREUR(208)
      endif
      lx(nx)=lw
      jx(nx)=jw
      x(nx)=rw
      if (ntx == 32) call DECFI(n2,32)
      end do
      call DECFI(n2,ntx)
      nn1=n1
      nn2=n2
      RETURN
      END subroutine TRIFI

!----------------------------------------------------------------------

!                          D E C F I

!----------------------------------------------------------------------
!     Decalage sur les arguments puis sur les n-uplets
!     des elements du tableau X(*) dans T(*,N) de numero N.
!     LONG est le nombre d'elements a decaler.
!----------------------------------------------------------------------
      subroutine DECFI(n,long)
      use PARTLIT
      implicit real*8 (a-h,o-z)

!----------------------------------------------------------------------
      if (long == 0) RETURN
      m=ntt(n)
      if (ntt(n) >= MXDIMT) call ERREUR(209)
      do i=0,long-1
      k=long-i
      nt=milt
      lw=lx(k)
      jw=jx(k)
      do j=idebt,16
         isix=jt(nt,n)
         if (jw > isix) then
            nt=nt+NPUIS2(j)
            l=1
         else
            if (jw == isix) then
            if (lw > lt(nt,n)) then
               nt=nt+NPUIS2(j)
               l=1
            else
               nt=nt-NPUIS2(j)
               l=0
            endif
            else
            nt=nt-NPUIS2(j)
            l=0
            endif
         endif
      end do
      nt=nt+l
      m=m+1
      jmax=m-nt
      if (jmax > 0) then
         do j=1,jmax
            ib=m-j
            ia=ib+k
            lt(ia,n)=lt(ib,n)
            t(ia,n)=t(ib,n)
            jt(ia,n)=jt(ib,n)
         end do
      endif
      if (jmax <  0) call ERREUR(210)
      lt(nt,n)=lw
      t(nt,n)=x(k)
      jt(nt,n)=jw
      m=nt
      end do
      ntt(n)=ntt(n)+long
      call RAZX
      RETURN
      END subroutine DECFI

!----------------------------------------------------------------------

!                          P U R G E

!----------------------------------------------------------------------
!     PURGE d'un tableau T(*,N) de numero N <>0 sur la precision PREC.
!----------------------------------------------------------------------
      subroutine PURGE(n)
      use PARTLIT
      implicit real*8 (a-h,o-z)

!----------------------------------------------------------------------
      if (n == 0) call ERREUR(221)
      call RAZ(n)
      jmax=ntr
      k=0
      do 1 j=1,jmax
      if (abs(r(j)) <  prec) goto 1
      k=k+1
      t(k,n)=r(j)
      lt(k,n)=lr(j)
      jt(k,n)=jr(j)
1     continue
      ntt(n)=k
      trig(n)=trigr
      RETURN
      END subroutine PURGE

!----------------------------------------------------------------------

!                          M U L S E R

!----------------------------------------------------------------------
!     Multiplication des tableaux d'indices N1, N2 rangee dans N3.
!     Si N3=0 le resultat n'est pas purge et reste dans R(*).
!     Si NATT=0, R(*) est remis a zero avant multiplication; si NATT=1,
!     R(*) n'est pas remis a zero (attente). FACTR est le facteur
!     multiplicatif du produit.
!     T(*,N3) ou R(*) = NATT*R(*) + FACTR*T(*,N1)*T(*,N2)
!     La limitation du produit est operee sur l'ordre des termes (MAXOP)
!     R(*) ou t(*,N3) est range par ordre croissant des n-uplets puis des
!     arguments.
!----------------------------------------------------------------------
      subroutine MULSER(nn1,nn2,nn3,natt,factr)
      use PARTLIT
      implicit real*8 (a-h,o-z)


      logical*1                           ::  trigw,tests,testd

!----------------------------------------------------------------------
      call CASH(nn1,.TRUE.,n1,nn2)
      if (nn1 == nn2) then
      n2=n1
      else
      call CASH(nn2,.TRUE.,n2,n1)
      end if
      trigw=trig(n1).eqv.trig(n2)
      if (natt == 0) then
      call RAZ(0)
      trigr=trigw
      else
      if (trigr.neqv.trigw) call ERREUR(222)
      endif
      tests=.FALSE.
      if (trigw.and.(.not.trig(n1))) tests=.TRUE.
      testd=.FALSE.
      if ((.not.trigw).and.trig(n1)) testd=.TRUE.
      call RAZX
      f=factr/2
      mx1=ntt(n1)
      mx2=ntt(n2)
      do 1 i=1,mx1
      iordre1=lt(i,n1)/idiv
      w=t(i,n1)*f
      do 2 j=1,mx2
         iordre2=lt(j,n2)/idiv
         iordre=iordre1+iordre2
         if (iordre <= maxop) then
            lww=lt(i,n1)+lt(j,n2)
            lw=lww
            ww=t(j,n2)*w
            rw=ww
            jw=jt(i,n1)+jt(j,n2)
            if (tests) rw=-ww
            call TRIARG
            jw=jt(i,n1)-jt(j,n2)
            rw=ww
            lw=lww
            if (testd) rw=-rw
            if (jw <  0) then
            jw=-jw
            if (.not.trigw) rw=-rw
            endif
            if (jw == 0.and..not.trigw) rw=0
            call TRIARG
         else
            if (j /=   1) goto 1
            goto 3
         endif
2     continue
1     continue
3     continue
      call DECAL(ntx)
      if (nn3 /=   0) then
      if (nn3 == nn1) n3=n1
      if (nn3 == nn2) n3=n2
      if (nn3 /=  nn1.and.nn3 /=  nn2) call CASH(nn3,.FALSE.,n3,0)
      call PURGE(n3)
      trig(n3)=trigr
      else
      n3=0
      endif
      n1=nn1
      n2=nn2
      n3=nn3
      RETURN
      END subroutine MULSER

!----------------------------------------------------------------------

!                          S O M S E R

!----------------------------------------------------------------------
!     Multiplication d'une serie par une constante (FACT) et addition
!     du resultat a l'accumulateur R(*).
!     Si N1 <>0: R(*) ou T(*,N2) = R(*) + FACT*T(*,N1).
!     Si N2 = 0 le resutat est stocke dans R(*), sinon il est range
!     dans T(*,N2). Si NATT=0, la table R(*) est remise a zero;
!     si NATT=1, R(*) n'est pas remise a zero (attente).
!     N1 peut etre nul, alors: R(*) ou T(*,N2) = FACT*R(*);
!     Dans ce cas R(*) n'est pas remise a zero (NATT inoperant).
!----------------------------------------------------------------------
      subroutine SOMSER(nn1,nn2,natt,fact)
      use PARTLIT
      implicit real*8 (a-h,o-z)

      logical*1                           ::  tg

!----------------------------------------------------------------------
      call CASH(nn1,.TRUE.,n1,0)
      call RAZX
      if (n1 == 0) then
      tg=trigr
      else
      tg=trig(n1)
      endif
      if (natt == 0.and.n1 /=   0) trigr=trig(n1)
      if (tg.neqv.trigr) call ERREUR(223)
      if (n1 == 0) then
      do i=1,ntr
         r(i)=fact*r(i)
      end do
      else
      if (natt == 0) then
         call RAZ(0)
      end if
      do i=1,ntt(n1)
         rw=fact*t(i,n1)
         lw=lt(i,n1)
         jw=jt(i,n1)
         call TRIARG
      end do
      call DECAL(ntx)
      end if
      if (nn2 /=   0) then
      if (nn2 == nn1) then
         n2=n1
      else
         call CASH(nn2,.FALSE.,n2,0)
      end if
      call PURGE(n2)
      trig(n2)=tg
      else
      n2=0
      endif
      nn1=n1
      nn2=n2
      RETURN
      END subroutine SOMSER

!----------------------------------------------------------------------

!                          S O M C T E

!----------------------------------------------------------------------
!     Addition d'une constante (FACT) et d'une serie a l'accumulateur $:
!     Si N1 <>0: R(*) ou T(*,N2) = R(*) + FACT + T(*,N1).
!     Si N2 = 0 le resultat est stocke dans R(*), sinon il est range
!     dans T(*,N2). Si NATT=0, la table R(*) est remise a zero; si NATT=1
!     R(*) n'est pas remise a zero (attente).
!     N1 peut etre nul; alors R(*) ou T(*,N2) = FACT + R(*)
!     Dans ce cas R(*) n'est pas remise a zero (NATT inoperant).
!----------------------------------------------------------------------
      subroutine SOMCTE(nn1,nn2,natt,fact)
      use PARTLIT
      implicit real*8 (a-h,o-z)

      logical*1                           :: tg

!----------------------------------------------------------------------
      call CASH(nn1,.TRUE.,n1,0)
      call RAZX
      if (natt == 0) then
      call RAZ(0)
      if (n1 /=   0) then
         trigr=trig(n1)
      else
         trigr=.TRUE.
      end if
      end if
      if (n1 == 0) then
      tg=trigr
      else
      tg=trig(n1)
      endif
      if (.not.tg) call ERREUR(224)
      if (n1 /=   0) then
      do i=1,ntt(n1)
         rw=t(i,n1)
         lw=lt(i,n1)
         jw=jt(i,n1)
         call TRIARG
      end do
      end if
      rw=fact
      lw=0
      jw=0
      call TRIARG
      call DECAL(ntx)
      if (nn2 /=   0) then
      if (nn1 == nn2) then
         n2=n1
      else
         call CASH(nn2,.FALSE.,n2,0)
      end if
      call PURGE(n2)
      trig(n2)=tg
      else
      n2=0
      endif
      nn1=n1
      nn2=n2
      RETURN
      END subroutine SOMCTE

!----------------------------------------------------------------------

!                          T R A N S

!----------------------------------------------------------------------
!     Tranfert de la serie T(*,N1) ou R(*) dans T(*,N2) ou R(*).
!----------------------------------------------------------------------
      subroutine TRANS(nn1,nn2)
      use PARTLIT
      implicit real*8 (a-h,o-z)

!----------------------------------------------------------------------
      if (nn1 == nn2) RETURN
      call CASH(nn1,.TRUE.,n1,0)
      call CASH(nn2,.FALSE.,n2,0)
      if (n1 == 0) then
      do i=1,ntr
         t(i,n2)=r(i)
         jt(i,n2)=jr(i)
         lt(i,n2)=lr(i)
      end do
      ntt(n2)=ntr
      trig(n2)=trigr
      else
      if(n2 == 0) then
         call RAZ(0)
         do i=1,ntt(n1)
            r(i)=t(i,n1)
            jr(i)=jt(i,n1)
            lr(i)=lt(i,n1)
         end do
         ntr=ntt(n1)
         trigr=trig(n1)
      else
         do i=1,ntt(n1)
            t(i,n2)=t(i,n1)
            jt(i,n2)=jt(i,n1)
            lt(i,n2)=lt(i,n1)
         end do
         ntt(n2)=ntt(n1)
         trig(n2)=trig(n1)
      end if
      end if
      nn1=n1
      nn2=n2
      RETURN
      END subroutine TRANS

!----------------------------------------------------------------------

!                          D E R I V E

!----------------------------------------------------------------------
!     Derivation par rapport au temps d'une serie N1.
!     Resultat dans N2
!----------------------------------------------------------------------
      subroutine DERIVE(nn1,nn2)
      use PARTLIT
      implicit real*8 (a-h,o-z)

      logical*1                           :: tg

      integer, dimension(NARG)            :: ifi
!----------------------------------------------------------------------

      call CASH(nn1,.TRUE.,n1,0)

      if (n1 == 0) call ERREUR(240)
      tg=trig(n1)

      imax=ntt(n1)
      do 1 i=1,imax
      iarg=jt(i,n1)
      call DECOMP(iarg,ifi)
      fact=0
         do ng=1,NARG
            fact=fact+ifi(ng)*freq(ng)
         end do
         if (tg) then
            r(i)=-fact*t(i,n1)
         else
            r(i)=fact*t(i,n1)
         end if
         jr(i)=jt(i,n1)
         lr(i)=lt(i,n1)
1     continue
      ntr=ntt(n1)
      trigr=.NOT.tg

      if (nn2 /= 0) then
         call CASH(nn2,.FALSE.,n2,0)
      call PURGE(n2)
      else
      n2=0
      endif

      nn1=n1
      nn2=n2
      END subroutine DERIVE

!----------------------------------------------------------------------

!                          P R I M I T

!----------------------------------------------------------------------
!     Primitive (integration par rapport au temps) d'une serie N1.
!     Resultat dans N2
!----------------------------------------------------------------------
      subroutine PRIMIT(nn1,nn2)
      use PARTLIT
      implicit real*8 (a-h,o-z)

      logical*1                           :: tg

      integer, dimension(NARG)            :: ifi
!----------------------------------------------------------------------

      call CASH(nn1,.TRUE.,n1,0)

      if (n1 == 0) call ERREUR(241)
      tg=trig(n1)

      imax=ntt(n1)
      do 1 i=1,imax
      iarg=jt(i,n1)
      call DECOMP(iarg,ifi)
      fact=0
         do ng=1,NARG
            fact=fact+ifi(ng)*freq(ng)
         end do
         if (fact == 0) call erreur (230)
         if (tg) then
            r(i)=t(i,n1)/fact
         else
            r(i)=-t(i,n1)/fact
         end if
         jr(i)=jt(i,n1)
         lr(i)=lt(i,n1)
1     continue
      ntr=ntt(n1)
      trigr=.NOT.tg

      if (nn2 /= 0) then
         call CASH(nn2,.FALSE.,n2,0)
      call PURGE(n2)
      else
      n2=0
      endif

      nn1=n1
      nn2=n2
      END subroutine PRIMIT

!----------------------------------------------------------------------

!                          D E R A R G

!----------------------------------------------------------------------
!     Derivee partielle par rapport a 1 argument de rang IPF
!     (fixe par #RANGARG) d'une serie de numero N1. Resultat dans N2.
!----------------------------------------------------------------------
      subroutine DERARG(nn1,nn2)
      use PARTLIT
      implicit real*8 (a-h,o-z)

      logical*1                           :: tg

      integer, dimension(NARG)            :: ifi
!----------------------------------------------------------------------

      call CASH(nn1,.TRUE.,n1,0)

      if (n1 == 0) call ERREUR(242)
      tg=trig(n1)

      imax=ntt(n1)
      do 1 i=1,imax
      iarg=jt(i,n1)
      call DECOMP(iarg,ifi)
      k=ifi(ipf)
         if (tg) then
            r(i)=-k*t(i,n1)
         else
            r(i)=k*t(i,n1)
         end if
         jr(i)=jt(i,n1)
         lr(i)=lt(i,n1)
1     continue
      ntr=ntt(n1)
      trigr=.NOT.tg

      if (nn2 /= 0) then
         call CASH(nn2,.FALSE.,n2,0)
      call PURGE(n2)
      else
      n2=0
      endif

      nn1=n1
      nn2=n2
      END subroutine DERARG

!----------------------------------------------------------------------

!                          D E R L I T

!----------------------------------------------------------------------
!     Derivee partielle par rapport a 1 parametre litteral de rang IPL
!     (fixe par #RANGLIT) d'une serie de numero N1. Resultat dans N2.
!----------------------------------------------------------------------
      subroutine DERLIT(nn1,nn2)
      use PARTLIT
      implicit real*8 (a-h,o-z)

      logical*1                           :: tg

      integer, dimension(NARG)            :: ifi
      integer, dimension(0:NLIT)          :: ilit
!----------------------------------------------------------------------

      call CASH(nn1,.TRUE.,n1,0)

      if (n1 == 0) call ERREUR(243)
      tg=trig(n1)

      call RAZX
      call RAZ(0)
      imax=ntt(n1)
      do 1 i=1,imax
      itas=lt(i,n1)
      call DETASSE(itas,ilit)
      k=ilit(ipl)
      if (k == 0) CYCLE
         jw=jt(i,n1)
         rw=k*t(i,n1)
         ilit(ipl)=k-1
         call TASSE(ilit,lw)
         call TRIARG
1     continue
      call DECAL(ntx)
      trigr=tg

      if (nn2 /= 0) then
         call CASH(nn2,.FALSE.,n2,0)
      call PURGE(n2)
      else
      n2=0
      endif

      nn1=n1
      nn2=n2

      END subroutine DERLIT

!----------------------------------------------------------------------

!                          S U B T I M E

!----------------------------------------------------------------------
!     Substitution du temps (TPS) et des valeurs numeriques des parametres
!     litteraux dans la serie de numero N1. Resultat dans RESUL
!----------------------------------------------------------------------
      subroutine SUBTIME(nn1,tps,resul)
      use PARTLIT
      implicit real*8      (a-h,o-z)

      logical*1                           :: tg

      integer, dimension(NARG)            :: ifi
      integer, dimension(0:NLIT)          :: ilit
!-----------------------------------------------------------------------
      call CASH(nn1,.TRUE.,n1,0)

      if (nn1 == 0) call ERREUR(244)
      tg=trig(n1)

      resul=0
      time=1
      imax=ntt(n1)
      resul=0
      do 1 i=1,imax

         itas=lt(i,n1)
      call DETASSE(itas,ilit)
         w=t(i,n1)
         do ng=1,NLIT
            if (ilit(ng) == 0) CYCLE
            w=w*VAR(ng)**ilit(ng)
         end do

         iarg=jt(i,n1)
      call DECOMP(iarg,ifi)
         fq=0
         phi=0
         do ng=1,NARG
            fq=fq+ifi(ng)*freq(ng)
            phi=phi+ifi(ng)*phase(ng)
         end do
         arg=fq*tps+phi
         if (tg) then
            resul=resul+w*cos(arg)
         else
            resul=resul+w*sin(arg)
         end if
1     continue

      nn1=n1
      END subroutine SUBTIME


!***********************************************************************
!***********************************************************************

!     ------------------- I M P R E S S I O N ----------------------

!***********************************************************************
!***********************************************************************

!----------------------------------------------------------------------

!                          I M P R I M

!----------------------------------------------------------------------
!     Commande globale d'impression qui redirige les ipressions vers
!     IMPRFLOT ou IMPRFRAC suivant la valeur de SIFRAC (sortie sous forme
!     fractionnaire ou non).
!     SIARGU : .FALSE. :la sortie est dans l'ordre croissant des ordres
!     (parametres litteraux puis arguments) : mode normal
!     SIARGU = .TRUE. : la sortie est dans l'ordre croissant des arguments
!     (arguments puis parametres litteraux) : mode 'argument'. Le nom des
!     series en sortie est precede d'un "underscore" (_)
!-----------------------------------------------------------------------
      subroutine IMPRIM(nom,n)
      use PARTLIT
      character*20                        :: nom,nomw
!----------------------------------------------------------------------
      ns=n
      if (SIARGU) then
         nomw='_'//nom
      if (nom == '$') then
         call AJOUTN(nomw)
      ns=numser(nser)
      call TRIFI(0,ns)
      else
         if (.not.xser) call ERREUR(111)
         call AJOUTN(nomw)
         ns=numser(nser)
         call TRIFI(n,ns)
      end if
      if (SIFRAC) call IMPRFRAC(ns)
         if (.NOT.SIFRAC) call IMPRFLOT(ns)
         call KILLN(nomw)
      else
         if (SIFRAC) call IMPRFRAC(ns)
         if (.NOT.SIFRAC) call IMPRFLOT(ns)
      end if

      END subroutine IMPRIM

!----------------------------------------------------------------------

!                          I M P R F L O T

!----------------------------------------------------------------------
!     Impression d'une serie en sinus (ou cosinus), avec des
!     coefficients en virgule flottante.
!     Si N <>0 la table T(*,N) est imprimee; Si N=0 la table
!     R(*) est imprimee. maxout est l'ordre maximum pour l'impression
!----------------------------------------------------------------------
      subroutine IMPRFLOT(n)
      use PARTLIT
      implicit real*8 (a-h,o-z)

      logical*1                           :: tg

      character zw*20,lignetri*3

      integer, dimension(NARG)            :: ifi
      integer, dimension(0:NLIT)          :: ilit

!----------------------------------------------------------------------
      call CASH(n,.TRUE.,n,0)
      if (n == 0) then
      zw='$'
      tg=trigr
      ntot=ntr
      else
      call QUELNOM(n,zw)
      tg=trig(n)
      ntot=ntt(n)
      end if
      write(*,*)
      write(*,'(A,A,A,i4)') ' Serie:   ',zw,' Nombre de termes: ',ntot
      if (tg) then
      lignetri='cos'
      else
      lignetri='sin'
      end if
      write(*,titre) lignetri
      if (n /=   0) then
         if (sishort) then
            imax=MIN(MAXIMPR,ntt(n))
         else
            imax=ntt(n)
         end if
      do 1 i=1,imax
         iarg=jt(i,n)
         itas=lt(i,n)
         call DETASSE(itas,ilit)
         if (ilit(0) > maxout) goto 1
         call DECOMP(iarg,ifi)
         write(*,zimp) i,ilit,t(i,n),ifi
1     continue
      RETURN
      else
         if (sishort) then
            imax=MIN(MAXIMPR,ntr)
         else
            imax=ntr
         end if
      do 2 i=1,imax
         iarg=jr(i)
         itas=lr(i)
         call DETASSE(itas,ilit)
         if (ilit(0) > maxout) goto 2
         call DECOMP(iarg,ifi)
         write(*,zimp) i,ilit,r(i),ifi
2     continue
      RETURN
      endif
      END subroutine IMPRFLOT

!----------------------------------------------------------------------

!                          I M P R F R A C

!----------------------------------------------------------------------
!     Impression d'une serie en sinus (ou cosinus), avec des coefficients
!     sous forme fractionnaire.
!     Si N <>0 la table T(*,N) est imprimee;
!     Si N=0 la table R(*) est imprimee. maxout est l'ordre maximum
!     pour l'impression.
!----------------------------------------------------------------------
      subroutine IMPRFRAC(n)
      use PARTLIT
      implicit real*8 (a-h,o-z)

      logical*1                           :: fail,tg

      character zw*20,lignetri*3

      integer, dimension(NARG)            :: ifi
      integer, dimension(0:NLIT)          :: ilit

!----------------------------------------------------------------------
      call CASH(n,.TRUE.,n,0)
      if (n == 0) then
      zw='$'
      tg=trigr
      ntot=ntr
      else
      call QUELNOM(n,zw)
      tg=trig(n)
      ntot=ntt(n)
      end if
      write(*,*)
      write(*,'(A,A,A,i4)') ' Serie:   ',zw,' Nombre de termes: ',ntot
      if (tg) then
      lignetri='cos'
      else
      lignetri='sin'
      end if
      write(*,titre) lignetri
      if (n /=   0) then
         if (sishort) then
            imax=MIN(MAXIMPR,ntt(n))
         else
            imax=ntt(n)
         end if
      do 1 i=1,imax
         iarg=jt(i,n)
         itas=lt(i,n)
         call DETASSE(itas,ilit)
         if (ilit(0) > maxout) goto 1
         call DECOMP(iarg,ifi)
         call FRACONT(t(i,n),inum,iden,fail)
         if (fail) then
            write(*,zimp1) i,ilit,t(i,n),ifi
         else
            write(*,zimp2) i,ilit,inum,iden,ifi
         end if
1     continue
      RETURN
      else
         if (sishort) then
            imax=MIN(MAXIMPR,ntr)
         else
            imax=ntr
         end if
      do 2 i=1,imax
         iarg=jr(i)
         itas=lr(i)
         call DETASSE(itas,ilit)
         if (ilit(0) > maxout) goto 2
         call DECOMP(iarg,ifi)
         call FRACONT(r(i),inum,iden,fail)
         if (fail) then
            write(*,zimp1) i,ilit,r(i),ifi
         else
            write(*,zimp2) i,ilit,inum,iden,ifi
         end if
2     continue
      RETURN
      endif
      END subroutine IMPRFRAC

!----------------------------------------------------------------------

!                          F R A C O N T

!----------------------------------------------------------------------
!     Calcul de la reduite d'un nombre F donne en entree;
!     en sortie: numerateur (INUM), denominateur (IDEN).
!     FAIL=.TRUE. si l'operation echoue (F est alors maintenu en REAL*8)
!     Les tolerances et les bornes du calcul sont initialisees au
!     1er passage dans la routine.
!----------------------------------------------------------------------
      subroutine FRACONT(f,inum,iden,fail)
      implicit real*8 (a-h,o-z)

      integer, parameter      :: MAXINT=2147483647
      logical*1               :: testol=.TRUE.,fail

!----------------------------------------------------------------------
      if (testol) then
      amx=MAXINT
      tol=1.d-14
      xtol=1/sqrt(100*tol)
      ytol=10*xtol
      testol=.FALSE.
      end if
      x=f
      if (x <  0) then
      is=-1
      x=-x
      else
      is=1
      endif
      fail=.FALSE.
      p0=1
      q0=0
      a=dint(x)
      p=a
      q=1
      b=x-a
      do 1 i=1,100
      if (b == 0) goto 2
      w=1/b
      a=dint(w)
      pp=p
      qq=q
      p=a*p+p0
      q=a*q+q0
      p0=pp
      q0=qq
      b=w-a
         if (abs(b) <= tol.or.abs(q) > xtol) goto 2
1     continue
      call ERREUR(225)
2     continue
      if (abs(q) > xtol) then
      p=p0
      q=q0
      endif
      a=abs(p/q-x)
      if (a > ytol) goto 3
      if (p > amx.or.q > amx) goto 3
      inum=p*is
      iden=q
      if (inum==0) goto 3           ! NEW 19-08-03
      RETURN
3     fail=.TRUE.
      RETURN
      END subroutine FRACONT

!***********************************************************************
!***********************************************************************

!     -- S E R I E S   D U  P R O B L E M E   D E S   2 - C O R P S --

!***********************************************************************
!***********************************************************************

!----------------------------------------------------------------------

!                          I N I H A N S

!----------------------------------------------------------------------
      subroutine INIHANS
      implicit real*8 (a-h,o-z)

      integer, parameter                  :: maxe=20,mx=maxe/2
      integer                             :: ordre,carac,smax,s,s1,s2
      real*8, dimension(0:mx)             :: f,p1,p2,p,sigma
      real*8, dimension(0:*)              :: h
      real*8, dimension(0:mx,0:maxe)      :: beta

!----------------------------------------------------------------------
      call DEVBET(mx,f)
      do j=1,mx
      beta(j,0)=0
      end do
      beta(0,0)=1
      do i=1,maxe
      kx=(maxe-i)/2
      call PROD(beta(0,i-1),f,kx,beta(0,i))
      end do
      RETURN

!----------------------------------------------------------------------

!                          H A N S E N

!----------------------------------------------------------------------
!     Calcul des coefficients de Hansen a l'ordre ORDRE en excentricite
!     (r/a)**n*exp i mv = Sigma(k) [X(n,m,k)*exp i kM]
!     Le calcul est realise a partir du formulaire de V. A. Brumberg.
!     Le coefficient X(n,m,k) de exp i kM est developpe dans H(*).
!     La caracteristique du terme est abs(m-k), c'est le degre
!     minimum en excentricite du developpement H(*).
!     L'initialisation des tables utiles a HANSEN est effectuee dans
!     le MAIN avec un CALL INIHANS. L'ordre maximum en excentricite est
!     fixe par MAXE. (ORDRE <= MAXE)
!----------------------------------------------------------------------
      entry HANSEN(N,M,K,ORDRE,H)

      mmk=m-k
      carac=abs(mmk)
      if (carac > ordre) call ERREUR(226)
      smax=(ordre-carac)/2
      if (k == 0) then
      am=carac
      a=(n+am+2)/2
      b=(n+am+3)/2
      c=1+am
      w=(-0.5)**carac*FACT(n+2,carac)/FACT(1,carac)
      call HYPERG(a,b,c,smax,f)
      call RAC(n+1.5d0,smax,p)
      call PROD(f,p,smax,sigma)
      do s=0,smax
         h(s)=w*sigma(s)
      end do
      else
      call DEVEL(-n-1,smax,f)
      imax=smax
      do s=0,smax
         sigma(s)=0
      end do
      do s=0,smax
         s1=s+max(0,-mmk)
         s2=s+max(0,mmk)
         call COEF_P(n,s1,m,k,imax,p1)
         call COEF_P(n,s2,-m,-k,imax,p2)
         call PROD(p1,p2,imax,p)
         do i=0,imax
            sigma(s+i)=sigma(s+i)+p(i)
         end do
         imax=imax-1
      end do
      call PROD(f,sigma,smax,p)
      do s=0,smax
         h(s)=0
      end do
      do s=0,smax
         do i=s,smax
            h(i)=h(i)+p(s)*beta(i-s,carac+2*s)
         end do
      end do
      end if
      RETURN
      END subroutine INIHANS

!----------------------------------------------------------------------

!                          D E V E L

!----------------------------------------------------------------------
!     Calcul des coefficients de F(*) = (1+x)**NA jusqu'au degre
!     x**IMAX (x est ici un parametre formel).
!----------------------------------------------------------------------
      subroutine DEVEL(na,imax,f)
      implicit real*8 (a-h,o-z)

      real*8, dimension(0:*)              :: f
!----------------------------------------------------------------------

      alpha=na
      f(0)=1
      do i=1,imax
      f(i)=(alpha-i+1)/i*f(i-1)
      end do
      RETURN
      END subroutine DEVEL

!----------------------------------------------------------------------

!                          R A C

!----------------------------------------------------------------------
!     Calcul des coefficients de F(*) = (1-x)**A jusqu'au degre
!     x**IMAX (x est ici un paramere formel).
!----------------------------------------------------------------------
      subroutine RAC(a,imax,f)
      implicit real*8 (a-h,o-z)

      real*8, dimension(0:*)              :: f

!----------------------------------------------------------------------
      f(0)=1
      do i=1,imax
      f(i)=-(a-i+1)/i*f(i-1)
      end do
      RETURN
      END subroutine RAC

!----------------------------------------------------------------------

!                          H Y P E R G

!----------------------------------------------------------------------
!     Calcul de la fonction hypergeometrique jusqu'a l'ordre IMAX;
!     le resultat est dans F(*).
!----------------------------------------------------------------------
      subroutine HYPERG(a,b,c,imax,f)
      implicit real*8 (a-h,o-z)

      real*8, dimension(0:*)              :: f

!----------------------------------------------------------------------
      if (c == 0) call ERREUR(227)
      f(0)=1
      do i=1,imax
      f(i)=(a+i-1)*(b+i-1)/((c+i-1)*i)*f(i-1)
      end do
      RETURN
      END subroutine HYPERG

!----------------------------------------------------------------------

!                          D E V B E T

!----------------------------------------------------------------------
!     Calcul des coefficients de F(*) = (1-rac(1-x**2))/x**2
!     jusqu'au degre 2*IMAX en x (x est ici un parametre formel).
!----------------------------------------------------------------------
      subroutine DEVBET(imax,f)
      implicit real*8 (a-h,o-z)

      real*8, dimension(0:*)              :: f

!----------------------------------------------------------------------
      alpha=0.5
      f(0)=0.5
      do i=1,imax
      f(i)=-(alpha-i)/(i+1)*f(i-1)
      end do
      RETURN
      END subroutine DEVBET

!----------------------------------------------------------------------

!                          F A C T

!----------------------------------------------------------------------
!     Calcul de  k(k+1)....(k+n-1) [n termes]
!----------------------------------------------------------------------
      function FACT(k,n)
      implicit real*8 (a-h,o-z)

!----------------------------------------------------------------------
      FACT=1
      do i=1,n
      FACT=FACT*(k+i-1)
      end do
      RETURN
      END function FACT

!----------------------------------------------------------------------

!                          C O E F _ P

!----------------------------------------------------------------------
!     Routine appelee par HANSEN.
!----------------------------------------------------------------------
      subroutine COEF_P(n,s,m,nu,imax,p)
      implicit real*8 (a-h,o-z)

      integer, parameter                  :: maxe=20,mx=maxe/2
      integer                             :: s,rstar,r
      real*8, dimension(0:mx)             :: f
      real*8, dimension(0:*)              :: p

!----------------------------------------------------------------------
      do i=0,imax
      p(i)=0
      end do
      u=nu
      ia=n-m+1
      ib=-ia
      if (ia >= 0) then
      rstar=min(s,ia)
      else
      rstar=s
      end if
      do r=0,rstar
      if (r == 0) then
         w=u**s/FACT(1,s)
      else
         w=((w*(ib+r-1))/r*(s-r+1))/nu
      end if
      call DEVEL(r-s,imax,f)
      do i=0,imax
         p(i)=p(i)+w*f(i)
      end do
      end do
      RETURN
      END subroutine COEF_P

!----------------------------------------------------------------------

!                          P R O D

!----------------------------------------------------------------------
!     Produit de deux polynomes P1(*)*P2(*) a l'ordre IMAX;
!     resultat dans P(*).
!----------------------------------------------------------------------
      subroutine PROD(p1,p2,imax,p)
      implicit real*8 (a-h,o-z)

      real*8, dimension(0:*)              :: p1,p2,p

!----------------------------------------------------------------------
      do i=0,imax
      p(i)=0
      end do
      do i=0,imax
      do j=i,imax
         p(j)=p(j)+p1(i)*p2(j-i)
      end do
      end do
      RETURN
      END subroutine PROD

!----------------------------------------------------------------------

!                          D E U X C O R P S

!----------------------------------------------------------------------
!     Evaluation des series du problemes des 2-corps a l'ordre ORDRE
!     en excentricite. La serie de numero ISER est chargee du
!     developpement litterale de (r/a)**n*cos(m*v) (XCOS=.TRUE.)
!     ou (r/a)**n*sin(m*v) (XCOS=.FALSE.). IPL est le rang de la
!     variable excentricite dans ILIT(*), tableau des parametres
!     litteraux, et IPF le rang de l'argument anomalie moyenne
!     dans IFI(*), tableau des arguments litteraux.
!----------------------------------------------------------------------
      subroutine DEUXCORPS(n,m,xcos,ordre,iser)
      use PARTLIT
      implicit real*8 (a-h,o-z)

      integer, parameter                  :: maxe=20,mx=maxe/2

      logical*1                           :: xcos,test

      integer                             :: ordre,carac

      integer, dimension(NARG)            :: ifi
      integer, dimension(0:NLIT)          :: ilit

      real*8, dimension(0:mx)             :: h

!----------------------------------------------------------------------
      call RAZX
      call RAZ(0)
      trigr=xcos
      do  i=1,NARG
      ifi(i)=0
      end do
      do  i=1,NLIT
      ilit(i)=0
      end do
      kmax=ordre+abs(m)
      if (xcos) then
      kmin=0
      else
      kmin=1
      end if
      do k=kmin,kmax
      carac=abs(k-m)
      if (carac <= ordre) call HANSEN(n,m,k,ordre,h)
      test=.TRUE.
10    continue
      jmax=ordre-carac
      do j=0,jmax,2
         ilit(ipl)=j+carac
         ifi(ipf)=k
         rw=h(j/2)
         if (.not.xcos.and..not.test) rw=-rw
         call COMP(ifi,jw)
         call TASSE(ilit,lw)
         call TRIARG
      end do
      if (test.and.k /=   0) then
         test=.not.test
         carac=abs(k+m)
         if (carac <= ordre) call HANSEN(n,-m,k,ordre,h)
         goto 10
      end if
      end do
      call DECAL(ntx)
      call CASH(iser,.FALSE.,iser,0)
      if (iser /=   0) call PURGE(iser)
      RETURN
      END subroutine DEUXCORPS

!***********************************************************************
!***********************************************************************

!  ------------ S A I S I E   D E S   S E R I E S ---------------

!***********************************************************************
!***********************************************************************

!----------------------------------------------------------------------

!                          S I N G L E

!----------------------------------------------------------------------
!     Subroutine d'une saisie simple (1 coefficient, 1 argument) de
!     la forme: A*x**n*sin (ou cos)(m*y) [x et y parametres formels].
!     Le resultat est dans T(*,ISER), (ou R(*) avec ISER=0).
!     XCOS=.TRUE.: serie en cosinus; XCOS=.FALSE.: serie en sinus.
!     IPL est le rang de la variable metrique [x] dans ILIT(*),
!     tableau des n-uplets litteraux, et IPF le rang de l'argument [y]
!     dans IFI(*), tableau des arguments litteraux.
!----------------------------------------------------------------------
      subroutine SINGLE(n,m,xcos,iser)
      use PARTLIT
      implicit real*8 (a-h,o-z)

      logical*1                           :: xcos

      integer, dimension(NARG)            :: ifi
      integer, dimension(0:NLIT)          :: ilit

!----------------------------------------------------------------------
      call RAZX
      call RAZ(0)
      trigr=xcos
      do  i=1,NARG
      ifi(i)=0
      end do
      do  i=1,NLIT
      ilit(i)=0
      end do
      ilit(ipl)=n
      ifi(ipf)=m
      rw=1
      call COMP(ifi,jw)
      call TASSE(ilit,lw)
      call TRIARG
      call DECAL(ntx)
      call CASH(iser,.FALSE.,iser,0)
      if (iser /=   0) call PURGE(iser)
      RETURN
      END subroutine SINGLE

!----------------------------------------------------------------------

!                          S A I S I E _ F

!----------------------------------------------------------------------
!     Subroutine de saisie d'une serie a partir d'un fichier.
!     Le nom du fichier de donnees qui contient la serie est NOMDATA.
!     Le resultat est dans T(*,ISER), (ou R(*) avec ISER=0).
!     XCOS=.TRUE.: serie en cosinus; XCOS=.FALSE.: serie en sinus.
!     Les donnees du fichier NOMDATA sont rangees dans l'ordre:
!     NLIT variable(s) metrique(s) chargee(s) dans ILIT(*),
!     tableau du n-uplet litteral;
!     1 coefficient associe au n-uplet (entier ou flottant);
!     NARG argument(s) charge(s) dans IFI(*), tableau des arguments
!     litteraux.
!----------------------------------------------------------------------
      subroutine SAISIE_F(nomdata,iser)
      use PARTLIT
      implicit real*8 (a-h,o-z)

      logical*1                           :: xcos

      character*80                        :: nomdata
      character*3                         :: trigono

      integer, dimension(NARG)            :: ifi
      integer, dimension(0:NLIT)          :: ilit


!----------------------------------------------------------------------
      OPEN (NLOG2,FILE=nomdata,STATUS='OLD',ERR=10)
      call RAZX
      call RAZ(0)
      read (NLOG2,'(A3)') trigono
      if (trigono=='COS') then
         xcos=.TRUE.
      else if (trigono=='SIN') then
         xcos=.FALSE.
      else
         print *,'Ligne trigonometrique non determinee'
      end if
      trigr=xcos
      do while (.TRUE.)
      read (NLOG2,FMT=*,END=20)
     &   (ilit(i), i=1,NLIT),rw,(ifi(i), i=1,NARG)
      call COMP(ifi,jw)
      call TASSE(ilit,lw)
      call TRIARG
      end do
20    continue
      CLOSE (NLOG2)
      call DECAL(ntx)
      call CASH(iser,.FALSE.,iser,0)
      if (iser /=   0) call PURGE(iser)
      RETURN
10    call ERREUR (128)
      END subroutine SAISIE_F

!----------------------------------------------------------------------

!                          S O R T I E _ F

!----------------------------------------------------------------------
!     Sortie d'une serie de numero N sur le fichier NOMDATA.
!     Le format de sortie (ZIMP3) pilote des enregistrements de series
!     utilisables ulterieurement a la lecture avec SAISIE_F
!----------------------------------------------------------------------

      subroutine SORTIE_F(nomdata,n)
      use PARTLIT
      implicit real*8 (a-h,o-z)

      character*80                        :: nomdata
      logical*1                           :: tg

      character lignetri*3

      integer, dimension(NARG)            :: ifi
      integer, dimension(0:NLIT)          :: ilit

!----------------------------------------------------------------------
      OPEN (NLOG2,FILE=nomdata)
      call CASH(n,.TRUE.,n,0)

      if (n == 0) then
      tg=trigr
      ntot=ntr
      else
      tg=trig(n)
      ntot=ntt(n)
      end if
      if (tg) then
      lignetri='COS'
      else
      lignetri='SIN'
      end if
      write (NLOG2,'(A3)') lignetri

      if (n /= 0) then
         imax=ntt(n)
      do 1 i=1,imax
         iarg=jt(i,n)
         itas=lt(i,n)
         call DETASSE(itas,ilit)
         if (ilit(0) > maxout) goto 1
         call DECOMP(iarg,ifi)
         write(NLOG2,zimp3) (ilit(k), k=1,NLIT),t(i,n),ifi
1     continue
      else
         imax=ntr
      do 2 i=1,imax
         iarg=jr(i)
         itas=lr(i)
         call DETASSE(itas,ilit)
         if (ilit(0) > maxout) goto 2
         call DECOMP(iarg,ifi)
         write(NLOG2,zimp3) (ilit(k), k=1,NLIT),r(i),ifi
2     continue
      endif
      close (NLOG2)
      END subroutine SORTIE_F

!----------------------------------------------------------------------

!                          O U T G R E G

!----------------------------------------------------------------------
!     Transformation d'une serie N1 en une serie de type 'GREGOIRE'
!     (Fourier et une seule ligne trigonometrique) apres
!     substitution des valeurs numeriques des parametres
!     litteraux dans la serie. Resultat dans N2
!----------------------------------------------------------------------
      subroutine OUTGREG(n1,n2)
      use PARTLIT
      implicit real*8      (a-h,o-z)

      logical*1                           :: tg

      integer, dimension(NARG)            :: ifi
      integer, dimension(0:NLIT)          :: ilit
!-----------------------------------------------------------------------
      if (n1 == 0) call ERREUR(247)
      call TRIFI(n1,n2)

      call RAZX
      call RAZ(0)
      imax=ntt(n2)

      do 1 i=1,imax
         itas=lt(i,n2)
      call DETASSE(itas,ilit)
         w=t(i,n2)
         do ng=1,NLIT
            if (ilit(ng) == 0) CYCLE
            w=w*VAR(ng)**ilit(ng)
         end do
         rw=w
         lw=0
         jw=jt(i,n2)
      call TRIARG
1     continue
      call DECAL(ntx)
      call PURGE(n2)

      END subroutine OUTGREG


!***********************************************************************
!***********************************************************************

!     ----  I N T E R P R E T E U R  D E   C O M M A N D E S  ----
!     ---- ( D E C O D A G E  D E S I N S T R U C T I O N S ) ----

!***********************************************************************
!***********************************************************************

!----------------------------------------------------------------------

!            I N T E R P R E T E

!-----------------------------------------------------------------------
!     Lecture d'une suite d'instructions dans ZREC. Les instructions
!     sont en sequence sur le fichier UNIT=NLOGY qui est lu jusqu'a
!     %e ou EOF.
!-----------------------------------------------------------------------
      subroutine INTERPRETE
      use PARTLIT
      implicit real*8      (a-h,o-z)

      logical              :: xboucle
      character*1          :: car
      character*2          :: first
      character*100        :: instruc
      character*101        :: zrec
!-----------------------------------------------------------------------
      xiter=.FALSE.
      xboucle=.TRUE.
10    continue
      read (NLOGY,FMT='(A)',END=20) zrec
!-----------------------------------------------------------------------
!     Elimination des commentaires entre crochets {...}
!-----------------------------------------------------------------------
      zrec(101:101)=' '
      zrecbis=zrec
      nd=INDEX(zrec,'{')
      do while (nd.ne.0)
         nf=INDEX(zrec,'}')
         if (nf == 0) call ERREUR(100)
         zrec=zrec(1:nd-1)//zrec(nf+1:)
         nd=INDEX(zrec,'{')
      end do
      if (INDEX(zrec,'}') /= 0) call ERREUR (100)
!-----------------------------------------------------------------------
!     Elimination des espaces.
!-----------------------------------------------------------------------
      nf=LEN_TRIM(zrec)
      if (nf == 0) goto 10
      nd=INDEX(zrec,' ')
      do while (nd < nf)
         zrec(nd:)=zrec(nd+1:)
         nf=LEN_TRIM(zrec)
         nd=INDEX(zrec,' ')
      end do
!-----------------------------------------------------------------------
!     Chargement des instructions dans la chaine INSTRUC a partir
!     de la suite ZREC. Le point-virgule (;) est le separateur
!     d'instructions. La chaine INSTRUC est chargee d'une seule
!     instruction. La derniere instruction de ZREC est optionnellement
!     suivie d'un point-virgule.
!     L'instruction est decodee par la routine ANAGAUCHE (membre
!     de gauche), qui appelle ANADROITE (membre de droite).
!     ZREC est exploree jusqu'a epuisement de la suite d'instructions.
!     (c'est a dire la purge de ZREC).
!-----------------------------------------------------------------------
      zrec(101:101)=';'
      nf=INDEX(zrec,';')
      do while (nf.ne.0)
         instruc=zrec(:nf-1)
!---------------------------------------------------------------------
!     Transformation des minuscules en majuscules sauf apres les commandes
!     %s (saisie d'une serie sur fichier) et %o (sortie de serie sur fichier)
!---------------------------------------------------------------------
         imax=nf-1
         nw=INDEX(instruc,'%s')+INDEX(instruc,'%S')
         nw=nw+INDEX(instruc,'%o')+INDEX(instruc,'%O')
         if (nw > 0) imax=nw+1
         do i=1,imax
            car=instruc(i:i)
            select case(car)
            case('a':'z')
               car=CHAR(ICHAR(car)-32)
               instruc(i:i)=car
            end select
         end do
!---------------------------------------------------------------------
         first=instruc(1:2)
         instrucbis=instruc
!-----------------------------------------------------------------------
!     Traitement de l'instruction dans la chaine: INSTRUC
!-----------------------------------------------------------------------
         if (instruc(1:1) /= ' ') call ANAGAUCHE(instruc)
         if (xfor) RETURN
         if (xiter) then
            if (first == '%B'.OR.first == '%b') then
               xboucle=.FALSE.
               ITER=ITER-1
               if (ITER == 0) then
                  xiter=.FALSE.
                  xboucle=.TRUE.
                  LIGNE=0
                  ITER=1
                  exit
               end if
               do i=0,LIGNE
                  BACKSPACE (NLOGY)
               end do
               exit
            end if
            if (xboucle) LIGNE=LIGNE+1
         end if
         zrec=zrec(nf+1:)
         nf=INDEX(zrec,';')
      end do
      goto 10
20    RETURN
      END subroutine INTERPRETE

!---------------------------------------------------------------------

!            A N A G A U C H E

!-----------------------------------------------------------------------
!     Decodage de la chaine INSTRUC contenant une seule instruction.
!     Si la routine renvoit au FORTRAN alors XFOR = .TRUE.
!     Sinon (cas standard) XFOR= .FALSE.
!-----------------------------------------------------------------------
      subroutine ANAGAUCHE(instruc)
      use PARTLIT
      implicit real*8      (a-h,o-z)

      logical              :: xwait,xposit
      character*1          :: car,cde,carcte,cwait
      character*2          :: first
      character*20         :: nom,nomc,nomw
      character*80         :: nomdata
      character*101        :: zrec
      character*(*)        :: instruc

      integer, dimension(0:NLIT)          :: ilit
      integer, dimension(NARG)            :: ifi
      character*20                        :: clist
      logical*1                           :: xrlit,xrarg

      integer, parameter              :: dimres=9
      character*20, dimension(dimres) :: nomres
      data nomres /'MAXOP','PREC','MAXOUT','SIFRAC','ITER',
     & 'MAXIMPR','RANGLIT','RANGARG','SIARGU'/

!-----------------------------------------------------------------------
      xfor=.FALSE.
      car=instruc(1:1)

      select case(car)
      case ('A':'Z','$','#')      ! Nom de serie ou cte
!-----------------------------------------------------------------------
!     Lecture d'un nom de serie ou de constante. Le nom est dans GAUCHE.
!     Le numero (eventuel) de la serie est dans NOP(0).
!     Si c'est une serie : XSERG = .TRUE.
!     Si c'est une constante : XCTEG = .TRUE. Le numero de la constante
!     est dans NCOP
!-----------------------------------------------------------------------
         if (car == '#') then
!---------------------------------------------------------------------
!     Recherche d'une constante reservee ou constante symbolique.
!     C'est une constante accessible par FORTRAN. Elle est precedee
!     du caractere (#)
!---------------------------------------------------------------------
            instruc=instruc(2:)
            call SYMBOLE(instruc,gauche,nx)
            do nres=1,dimres
               if (TRIM(gauche) == nomres(nres)) then
                  n=INDEX(instruc,'=')
                  if (n == 0) call ERREUR(110)
                  instruc=instruc(n+1:)
                  call NOMBRE(instruc,dx,nx)
                  ndx=IDINT(dx)
                  select case(nres)
                  case (1)
                     MAXOP=ndx
                     if (MAXOP > MAXORD) call ERREUR(250)
                     write (*,'(1X,A,I4)')
     &      'LA TRONCATURE ANALYTIQUE DES SERIES EST (MAXOP): ',MAXOP
                  case (2)
                     PREC=dx
                     write (*,'(1X,A,E10.2)')
     &      'LA TRONCATURE DE PURGE DES SERIES EST (PREC): ',PREC
                  case (3)
                     MAXOUT=ndx
                  case (4)
                     if (ndx < 0 .OR. ndx > 1) call ERREUR(109)
                     if (ndx==1) SIFRAC=.TRUE.
                     if (ndx==0) SIFRAC=.FALSE.
                  case (5)
                     if (ndx <= 0) call ERREUR(108)
                     ITER=ndx
                   case (6)
                     MAXIMPR=ndx
                   case (7)
                     if (ndx < 0 .OR. ndx > NLIT) call ERREUR(116)
                     ipl=ndx
                  case (8)
                     if (ndx < 0 .OR. ndx > NARG) call ERREUR(116)
                     ipf=ndx
                  case (9)
                     if (ndx < 0 .OR. ndx > 1) call ERREUR(109)
                     if (ndx==1) SIARGU=.TRUE.
                     if (ndx==0) SIARGU=.FALSE.
                  case default
                     call ERREUR (106)
                  end select

                  RETURN
               end if
            end do
            call ERREUR (106)
         end if

         if (car == '$') then
            nop(0)=0
            xserg=.TRUE.
            xcteg=.FALSE.
         else
            call SYMBOLE(instruc,gauche,nx)
!---------------------------------------------------------------------
!     Recherche d'une constante ou d'un nom de serie
!---------------------------------------------------------------------
            call TRIC(gauche)
            xcteg=xcte
            if (xcte) ncop=ncte
            call TRIN(gauche)
            xserg=xser
            if (xser) nop(0)=numser(nser)
         end if

         n=INDEX(instruc,'=')
         if (n == 0) call ERREUR(110)
         instruc=instruc(n+1:)

         if (instruc(1:1) /= '%') then
            call SICONST(instruc)
            if (xtot) then
!---------------------------------------------------------------------
!     Analyse du membre droit dans le cas exclusif ou il s'agit
!     d'operations sur des constantes
!     Toute cette sequence se poursuit jusqu'a: 25 continue
!---------------------------------------------------------------------
               if (xserg) goto 25
               carcte=instruc(1:1)
               xwait=.FALSE.
               s=0
               do while (carcte /= ' ')
                  xposit=.TRUE.
                  select case (carcte)
                  case('+','-')
                     if (carcte == '-') xposit=.FALSE.
                     instruc=instruc(2:)
                     carcte=instruc(1:1)
                  end select
                  select case (carcte)
                  case ('0':'9','.')
                     call NOMBRE(instruc,c,mr)
                  case ('A':'Z')
                     call SYMBOLE(instruc,nomc,mr)
                     call TRIC(nomc)
                     if (.NOT.xcte) call ERREUR(138)
                     c=cte(ncte)
                  case default
                     call ERREUR(138)
                  end select
                  if (.NOT.xposit) c=-c
                  carcte=instruc(mr:mr)
                  if (xwait) then
                     select case (cwait)
                     case ('*')
                        p=cw*c
                     case ('/')
                        if (ABS(c) <= 1.D-20) call erreur(139)
                        p=cw/c
                     end select
                  else
                     p=c
                  end if
                  select case (carcte)
                  case ('*','/')
                     cw=c
                     if (xwait) call ERREUR(170)
                     xwait=.TRUE.
                     cwait=carcte
                     p=0
                     instruc=instruc(mr+1:)
                     if (instruc == ' ') call ERREUR(170)
                  case default
                     xwait=.FALSE.
                     instruc=instruc(mr:)
                  end select
                  s=s+p
                  carcte=instruc(1:1)
               end do
               if (.NOT.xcteg) then
                  call AJOUTC(gauche,s)
               else
                  cte(ncop)=s
               end if
               RETURN
            end if
25          continue
            call ANADROITE(instruc)
         else
!---------------------------------------------------------------------
!     Derriere le signe (=) on trouve un operateur (%), sinon on
!     appelle ANADROITE pour l'analyse du membre droit.
!     Liste des operateurs admissibles:
!     s (S) : Saisie d'une serie
!     o (O) : Sortie d'une serie (format de saisie)
!     d (D) : Derivation d'une serie
!     p (P) : Primitive d'une serie
!     y (Y) : Transformation en une serie de type 'GREGOIRE'
!     z (Z) : Subtitution du temps
!     h (H) : Recherche dans une serie d'un monome ou d'un argument ou les deux
!---------------------------------------------------------------------
            cde=instruc(2:2)
            select case(cde)
            case ('O')
               if (xcteg) call ERREUR(131)
               if (xserg) then
                  nout=nop(0)
               else
                  call ERREUR(111)
               end if
               instruc=instruc(3:)
               nz=INDEX(instruc,' ')
               if (nz == 1) call ERREUR(129)
               nomdata=instruc(:nz-1)
               call SORTIE_F(nomdata,nout)
            case ('S')
!---------------------------------------------------------------------
!     Saisie d'une serie sur le fichier NOMDATA
!     (avec ou sans action sur les arguments avec le test de NEWARGU)
!---------------------------------------------------------------------
               if (xcteg) call ERREUR(131)
               if (xserg) then
                  nout=nop(0)
               else
                  call AJOUTN(gauche)
                  nout=numser(nser)
               end if
               instruc=instruc(3:)
               nz=INDEX(instruc,' ')
               if (nz == 1) call ERREUR(129)
               nomdata=instruc(:nz-1)
               call SAISIE_F(nomdata,nout)
!----------------------------------------------------------------------
!     Recherche du et/ou des termes d'une serie de parametres litteraux
!     et/ou d'arguments donnes.
!     Les termes sont affiches sur l'unite de sortie
!----------------------------------------------------------------------
            case ('H')
               if (xcteg) call ERREUR(131)
               if (xserg) then
                  nout=nop(0)
               else
                  call ERREUR(111)
               end if
               instruc=instruc(3:)
               nb=INDEX(instruc,'(')
               if (nb==0) call ERREUR (251)
               instruc=instruc(nb+1:)
               nb=INDEX(instruc,')')
               if (nb==0) call ERREUR (251)
               if (nb==1) then
                  xrlit=.FALSE.
               else
                  xrlit=.TRUE.
                  clist=instruc(1:nb-1)
                  read (clist,*,IOSTAT=istat) (ilit(k), k=1,NLIT)
                  if (istat/=0) then
                     if (istat>0) print *,'ERREUR DANS LA LISTE'
                     if (istat<0) print *,'LISTE INCOMPLETE'
                     call ERREUR(251)
                  end if
                  call tasse(ilit,lw)
               end if
               instruc=instruc(nb+1:)
               nb=INDEX(instruc,'(')
               if (nb==0) call ERREUR(251)
               instruc=instruc(nb+1:)
               nb=INDEX(instruc,')')
               if (nb==0) call ERREUR(251)
               if (nb==1) then
                  xrarg=.FALSE.
               else
                  xrarg=.TRUE.
                  clist=instruc(1:nb-1)
                  read (clist,*,IOSTAT=istat) (ifi(k), k=1,NARG)
                  if (istat/=0) then
                     if (istat>0) print *,'ERREUR DANS LA LISTE'
                     if (istat<0) print *,'LISTE INCOMPLETE'
                     call ERREUR(251)
                  end if
                  call COMP(ifi,jw)
               end if
               if (.NOT.xrlit .AND. .NOT.xrarg) call ERREUR(251)
               call QUELNOM(nout,nomw)
               write (*,'(1X,A)') 'Recherche dans la serie: '//nomw
               imax=ntt(nout)
               if (xrlit.AND.xrarg) then
                  do i=1,imax
                     if (jt(i,nout)==jw.AND.lt(i,nout)==lw) then
                        write (*,zimp) i,ilit,t(i,nout),ifi
                     end if
                  end do
               end if
               if (xrlit.AND..NOT.xrarg) then
                  do i=1,imax
                     if (lt(i,nout)==lw) then
                        iarg=jt(i,nout)
                        call DECOMP(iarg,ifi)
                        write (*,zimp) i,ilit,t(i,nout),ifi
                     end if
                  end do
               end if
               if (xrarg.AND..NOT.xrlit) then
                  do i=1,imax
                     if (jt(i,nout)==jw) then
                        itas=lt(i,nout)
                        call DETASSE(itas,ilit)
                        write (*,zimp) i,ilit,t(i,nout),ifi
                     end if
                  end do
               end if

!---------------------------------------------------------------------
!     Derivation d'une serie (et une seule)
!---------------------------------------------------------------------
            case ('D','P','Y')
               if (xcteg) call ERREUR(131)
               if (xserg) then
                  nout=nop(0)
               else
                  call AJOUTN(gauche)
                  nout=numser(nser)
               end if
               instruc=instruc(3:)
               nb=LEN_TRIM(instruc)
               call SYMBOLE(instruc,nom,nx)
               call TRIN(nom)
               if (nom == '$') then
                  nin=0
               else
                  nin=numser(nser)
               end if
               select case (cde)
               case ('D')
                  call DERIVE(nin,nout)
               case ('P')
                  call PRIMIT(nin,nout)
               case ('Y')
                  call OUTGREG(nin,nout)
               end select
            case ('Z')
               if (xserg) call ERREUR(119)
               instruc=instruc(3:)
               call SYMBOLE(instruc,nom,nx)
               call TRIN(nom)
               if (nom == '$') xser=.TRUE.
               if (.NOT.xser) call ERREUR(119)
               if (nom == '$') then
                  nin=0
               else
                  nin=numser(nser)
               end if
               if (instruc(nx:nx) /= ',') call ERREUR(119)
               instruc=instruc(nx+1:)
               call SYMBOLE(instruc,nom,nx)
               call TRIC(nom)
               if (.NOT.xcte) call ERREUR(119)
               tps=cte(ncte)
               call SUBTIME(nin,tps,s)
               if (.NOT.xcteg) then
                  call AJOUTC(gauche,s)
               else
                  cte(ncop)=s
               end if
            case default
               call ERREUR(113)
            end select
         end if
!-----------------------------------------------------------------------
!     Lecture d'une commande:
!     e(E) : fin des operations
!     f(F) : Retour au FORTRAN
!     m(M) : Redirection d'un message a la console
!         m Voici un joli message...
!     i(I) XXX, YYY ...Impression des series.
!         Impression des series de noms XXX, YYY, ...
!         jusqu'a epuisement de la liste. Controle de l'impression
!         avec MAXOUT
!     j(J) XXX, YYY ... dito i(I). Controle de l'impression avec MAXIMPR
!     t(T) XXX, YYY ... Elimination (tuer) les series XXX, YYY, ...
!        jusqu'a epuisement de la liste
!     x(X) : Iteration
!     b(B) : Fin de boucle
!     +    : Incrementation de l'INDICE N
!     g(G) : Appel d'une routine utilisateur
! ATTENTION ! Dans le cas d'une liste multiple, chaque element doit
! etre separe par une virgule (,)
!-----------------------------------------------------------------------
      case ('%')                                      ! Symbole de commande
         cde=instruc(2:2)
         select case (cde)
         case ('E')
            print *, 'FIN DES OPERATIONS'
            STOP
         case ('F')
            xfor=.TRUE.
            RETURN
         case ('B')
            if (.NOT.xiter) call ERREUR(136)
         case ('X')
!---------------------------------------------------------------------
!     La sequence qui suit est un garde fou pour tester la fin de
!     boucle
!---------------------------------------------------------------------
            lw=0
            read (NLOGY,FMT='(A)',END=30) zrec
            first=zrec(1:2)
            do while (first /= '%B'.AND.first /= '%b')
               lw=lw+1
               read (NLOGY,FMT='(A)',END=30) zrec
               first=zrec(1:2)
               if (first == '%X'.OR.first == '%x') goto 30
            end do
            do i=0,lw
               BACKSPACE (NLOGY)
            end do
!---------------------------------------------------------------------
            xiter=.TRUE.
            RETURN
30          call ERREUR(137)
!-----------------------------------------------------------------------
         case ('+')
            instruc=instruc(3:)
            nb=LEN_TRIM(instruc)
            if (nb == 0) then
               INDICE=INDICE+1
            else
               call NOMBRE(instruc,dx,nx)
               ndx=IDINT(dx)
               INDICE=INDICE+ndx
            end if
         case ('G')
            instruc=instruc(3:)
            call EXPLORE(instruc)
         case ('M')
            instruc=zrecbis(4:)
            print *,instruc
         case ('I','J','T')
            instruc=instruc(3:)
            nb=LEN_TRIM(instruc)
            if (nb == 0) call ERREUR(140)
            do while (nb.ne.0)
               call SYMBOLE(instruc,nom,nx)
               call TRIC(nom)
               call TRIN(nom)
               n1=numser(nser)
               select case (cde)
               case ('I')
                  sishort=.FALSE.
                  if (nom == '$') then
                     nw=0
                     call IMPRIM(nom,nw)
                  else
                     if (xser) then
                        call IMPRIM(nom,n1)
                     else
                        if (.NOT.xcte) call ERREUR(171)
                        print *,nom//'=',cte(ncte)
                     end if
                  end if
               case ('J')
                  sishort=.TRUE.
                  if (nom == '$') then
                     nw=0
                     call IMPRIM(nom,nw)
                  else
                     if (.not.xser) call ERREUR(111)
                     call IMPRIM(nom,n1)
                  end if
               case ('T')
                  if (.not.xcte.and..not.xser) call ERREUR(112)
                  if (xcte) call KILLC(nom)
                  if (xser) call KILLN(nom)
               case default
                  call ERREUR(113)
               end select
!     Recherche d'un separateur (,) dans le cas d'une liste
               ny=INDEX(instruc,',')
               if (ny < nx ) then
                  nb=INDEX(instruc(nx:),',')
               else
                  nb=ny
               end if
               instruc=instruc(nx+1:)
            end do
         case default
            call ERREUR(113)
         end select
      case default
         call ERREUR(114)
      end select                              ! Fin des recherches (nom ou %)
      END subroutine ANAGAUCHE

!-----------------------------------------------------------------------


!            A N A D R O I T E

!-----------------------------------------------------------------------
!     Analyse du membre de droite de l'instruction dans INSTRUC.
!     L'instruction est exploree de gauche a droite.
!     NCYCLE est l'indicateur du cycle d'operation;
!     NTYPE est l'indicateur du type operatoire (cf. CYCLOPE).
!     NCYCLE=NCYCLE+1 a chaque nouveau groupe d'operations.
!     NTYPE=NTYPE+1 avec une constante; NTYPE=NTYPE+2 avec une serie.
!     La valeur de la constante (eventuelle) dans chaque cycle est
!     chargee dans COP. Les numeros des series du groupe operatoire
!     sont charges dans NOP: NOP(0) serie du membre gauche, NOP(1:2)
!     serie(s) du membre droit. L'appel a CYCLOPE est effectue avec
!     XOUT=.FALSE. a la fin de chaque cycle, et avec XOUT=.TRUE. a la fin
!     de l'instruction. La soustraction de deux series est prise en
!     compte avec SPLUS=.FALSE., (changement du signe de COP)
!     sinon l'addition est effectuee sans changement de signe.
!-----------------------------------------------------------------------
      subroutine ANADROITE(instruc)
      use PARTLIT
      implicit real*8      (a-h,o-z)

      logical              :: ctecy,splus
      character*1          :: car
      character*20         :: droite
      character*(*)        :: instruc
!-----------------------------------------------------------------------
      un=1
      ncycle=1
      splus=.TRUE.
1     continue
      ntype=0
      is=1
      c=un
      ctecy=.FALSE.
      car=instruc(1:1)
      do itour=1,3
         select case(car)
         case ('$')
            if (ncycle > 1) call ERREUR(120)
            ntype=ntype+2
            nop(is)=0
            is=is+1
            mr=2
         case ('A':'Z')
            call SYMBOLE(instruc,droite,mr)
            call TRIC(droite)
            if (xcte) then
               if (ctecy) call ERREUR(121)
               c=cte(ncte)
               ctecy=.TRUE.
               ntype=ntype+1
            else
               call TRIN(droite)
               if (.not.xser) call ERREUR(122)
               ntype=ntype+2
               if (ntype > 5) call ERREUR(123)
               nop(is)=numser(nser)
               is=is+1
            end if
         case ('0':'9','+','-','.')
            if (ctecy) call ERREUR(121)
            call NOMBRE(instruc,c,mr)
            if (mr==2.AND.car=='-') then
               mr=mr-1
            end if
            ctecy=.TRUE.
            ntype=ntype+1
         case default
            call ERREUR(126)
         end select
         car=instruc(mr:mr)
         select case (car)
         case ('*')
         case ('+','-')
            if (splus) then
               cop=c
            else
               cop=-c
            endif
            if (car == '-') then
               splus=.FALSE.
            else
               splus=.TRUE.
            end if
            call CYCLOPE(ncycle,ntype,.FALSE.)
            ncycle=ncycle+1
            instruc=instruc(mr+1:)
            goto 1
         case (' ')
            if (splus) then
               cop=c
            else
               cop=-c
            endif
            call CYCLOPE(ncycle,ntype,.TRUE.)
            RETURN
         case default
            call ERREUR(124)
         end select
         instruc=instruc(mr+1:)
         car=instruc(1:1)
      end do
      call ERREUR(125)
      END subroutine ANADROITE

!-----------------------------------------------------------------------

!            C Y C L O P E

!-----------------------------------------------------------------------
!     Appel des routines d'addition et/ou de multiplication.
!     Si XOUT=.TRUE., le cycle est termine, sinon XOUT=.FALSE.
!     NCYCLE est le numero du cycle operatoire (cycle: groupe
!     d'operations successives A*B*C). Le type NTYPE definit les
!     operations. (Notations C=cte; S=serie)
!     NTYPE=1          C
!     NTYPE=2          S
!     NTYPE=3          C*S
!     NTYPE=4          S*S
!     NTYPE=5          C*S*S
!     Les resultats des operations des cycles successifs sont
!     accumules dans $.
!-----------------------------------------------------------------------
      subroutine CYCLOPE(ncycle,ntype,xout)
      use PARTLIT
      implicit real*8      (a-h,o-z)

      logical              :: xout
!-----------------------------------------------------------------------

      if (xout) then
         if (ntype == 1.and.ncycle == 1) then
            if (xcteg) then
               cte(ncop)=cop
            else
               if (xserg) call ERREUR(130)
               call AJOUTC(gauche,cop)
            endif
            RETURN
         else
            if (xcteg) call ERREUR(131)
            if (xserg) then
               nout=nop(0)
            else
               call AJOUTN(gauche)
               nout=numser(nser)
            endif
         endif
      else
         nout=0
      endif
      if (ncycle == 1) then
         select case (ntype)
         case (1)
            nw=0
            call SOMCTE(nw,nout,0,cop)
         case (2)
            call TRANS(nop(1),nout,.TRUE.)
         case (3)
            call SOMSER(nop(1),nout,0,cop)
         case (4:5)
            if (nop(1) == 0.or.nop(2) == 0) call ERREUR(132)
            call MULSER(nop(1),nop(2),nout,0,cop)
         case default
            call ERREUR (133)
         end select
      else
         select case (ntype)
         case (1)
            nw=0
            call SOMCTE(nw,nout,1,cop)
         case (2:3)
            call SOMSER(nop(1),nout,1,cop)
         case (4:5)
            if (nop(1) == 0.or.nop(2) == 0) call ERREUR(132)
            call MULSER(nop(1),nop(2),nout,1,cop)
         case default
            call ERREUR(133)
         end select
      end if
1     continue
      END subroutine CYCLOPE


!----------------------------------------------------------------------

!                          S Y M B O L E

!----------------------------------------------------------------------
!     Transformation d'une chaine de caracteres en un nom de serie
!     (majuscules, nombre, "_", ou "$"). Entree: CHAINE; Sortie: NOM
!     (chaine <=10 carac.), NX, rang du 1er caractere de chaine
!     n'appartenant pas au nom. Les minuscules sont transformees en
!     majuscules.
!----------------------------------------------------------------------
      subroutine SYMBOLE(chaine,nom,nx)
      use PARTLIT
      implicit real*8 (a-h,o-z)

      character*1                         :: car
      character*20                        :: nom
      character*(*)                       :: chaine

!----------------------------------------------------------------------
      l=len(chaine)
      if (l == 0) call ERREUR(140)
      do i=1,l
      car=chaine(i:i)
      select case(car)
      case('a':'z')
         car=char(ichar(car)-32)
         chaine(i:i)=car
      case ('A':'Z','_','$','1':'9')
      case default
         nx=i
         go to 100
      end select
      end do
      nx=l+1
100   if (nx > 21) call ERREUR(141)
      nom=chaine(:nx-1)
      RETURN
      END subroutine SYMBOLE

!----------------------------------------------------------------------

! [G]                         N O M B R E

!----------------------------------------------------------------------
!     Transformation d'une chaine de caracteres en un flottant
!     (double precision). Entree: CHAINE; Sortie: DX (flottant),
!     NX: rang du 1er caractere de chaine n'appartenant pas au flottant.
!----------------------------------------------------------------------
      subroutine NOMBRE(chaine,dx,nx)
      use PARTLIT
      implicit real*8 (a-h,o-z)

      logical*1                           :: point,expo,digit
      character*1                         :: car
      character*50                        :: cx
      character*(*)                       :: chaine

!----------------------------------------------------------------------
      point=.FALSE.
      expo=.FALSE.
      digit=.FALSE.
      j=0
      l=len(chaine)
      if (l == 0) call ERREUR(140)
      do i=1,l
      car=chaine(i:i)
      select case(car)
      case ('0':'9')
         digit=.TRUE.
      case ('.')
         if (point) call ERREUR (142)
         point=.TRUE.
      case ('d','D','e','E')           !!!!! a revoir
         if (expo) call ERREUR(142)
         expo=.TRUE.
         j=i
         digit=.FALSE.
      case ('+','-')
         if (i > 1) then
            if (.not.expo) then
            if (.not.digit) call ERREUR(142)
            nx=i
            goto 100
            else
            if ((i-j) > 1) then
               if (.not.digit) call ERREUR(142)
               nx=i
               goto 100
            endif
            endif
         endif
      case default
         if (i == 1) call ERREUR(142)
         nx=i
         go to 100
      end select
      end do
      nx=l+1
100   if (nx > 41) call ERREUR(143)
      cx=chaine(:nx-1)
      read (cx,'(F40.0)') dx
      RETURN
      END subroutine NOMBRE

!---------------------------------------------------------------------

!                             E X P L O R E

!---------------------------------------------------------------------
!     Routine d'exploration d'une ligne de commande (%G) qui engendre
!     un appel a une subroutine de nom connu.
!     La ligne de commande a la forme suivante:
!     (%G) SUB, PAR1, PAR2, PAR3 , ...
!     SUB est le nom d'une subroutine presente dans la data SUB. PAR_i
!     sont des parametres :noms de tableaux, noms de constantes symboliques,
!     constantes numeriques ou $.
!     La routine classe dans l'ordre de la liste la nature des parametres
!     dans PARAM(*), les numeros (eventuels) de series ou de constantes
!     dans NSTOC(*) et les valeurs (eventuelles) des constantes dans
!     CSTOC(*). Si la constante/serie est de nom inconnu PARAM=5.
!     Sa valeur est stockee dans NRESUL(*)
!     C'est le nom d'un symbole ou transporter le resultat de la subroutine.
!     S'il y a lieu, la constante/serie sera cree par un :
!     call AJOUTN(NRESUL(*)) pour une serie
!     call AJOUTC(NRESUL(*),C) pour une constante C.
!
!     Signification des parametres de PARAM:
!     0 : $
!     1 : serie existante
!     2 : constante numerique
!     3 : constante symbolique
!     4 : subroutine
!     5 : resultat de nature inconnu
!
!     Le rang dans SUB de la subroutine se trouve dans NSTOC(1). Ce rang
!     permet l'appel de la subroutine. Sa liste d'appel est constituee
!     avec les 3 tableaux PARAM, CSTOC et NSTOC pour son utilisation.
!-----------------------------------------------------------------------
      subroutine EXPLORE(instruc)
      use PARTLIT
      implicit real*8      (a-h,o-z)

      logical*1                        :: xcos
      character*1                      :: car
      character*20                     :: nom
      character*(*)                    :: instruc
      character*101                 :: zrec

      integer, parameter               :: dimstoc=5
      integer, dimension (dimstoc)     :: param,nstoc
      real*8, dimension(dimstoc)       :: cstoc
      character*20, dimension(dimstoc) :: nresul
      integer, parameter               :: nsub=9
      character*20, dimension(nsub)    :: sub
      data sub /'CHERCHEC','CONSER','HANSENC','HANSENS',
     & 'SINGLEC','SINGLES','DERARG','DERLIT','SINCOS'/
!-----------------------------------------------------------------------
!     Exploration de la ligne de commande
!-----------------------------------------------------------------------
      nb=LEN_TRIM(instruc)
      if (nb == 0) call ERREUR(140)
      ir=0
      do while (nb.ne.0.AND.ir < dimstoc)
         ir=ir+1
         if (ir == 1) then
            call SYMBOLE(instruc,nom,nx)
            do n=1,nsub
               if (TRIM(nom) == sub(n)) goto 20
            end do
            if (TRIM(NOM) == 'END') then
               CLOSE (NLOGY)
               NLOGY=NLOG1
               RETURN
            end if

            subutil=nom
            NLOGY=NLOG10
            OPEN (NLOG10,FILE=FILSUB)

100         continue
            read (NLOGY,FMT='(A)',END=200) zrec
            nw=INDEX(zrec,TRIM(nom))
            if (nw == 0) goto 100
            RETURN
200         call ERREUR(149)

20          param(ir)=4                   ! subroutine
            nstoc(ir)=n                   ! numero de la subroutine dans SUB
            goto 10
         end if
         car=instruc(1:1)
         select case (car)
         case ('0':'9','+','-','.')
            call NOMBRE(instruc,c,nx)
            param(ir)=2                   ! constante numerique
            cstoc(ir)=c                   ! valeur dans CSTOC
         case ('$')
            call SYMBOLE(instruc,nom,nx)
            param(ir)=0                   ! accumulateur $
            nstoc(ir)=0                   ! numero de la serie (ici $)
         case ('A':'Z')
            call SYMBOLE(instruc,nom,nx)
            call TRIC(nom)
            if (xcte) then
               param(ir)=3                ! constante symbolique
               cstoc(ir)=cte(ncte)        ! valeur dans CSTOC
               nstoc(ir)=ncte
               goto 10
            end if
            call TRIN(nom)
            if (xser) then
               param(ir)=1                ! serie
               nstoc(ir)=numser(nser)
               goto 10
            else
               nresul(ir)=nom
               param(ir)=5                ! Resultat de nature inconnue
            end if
         case default
            call ERREUR(145)
         end select
!     Recherche d'un separateur (,) dans le cas d'une liste
10       continue
         ny=INDEX(instruc,',')
         if (ny < nx ) then
            nb=INDEX(instruc(nx:),',')
         else
            nb=ny
         end if
         instruc=instruc(nx+1:)
      end do
!-----------------------------------------------------------------------
!     Appel de la subroutine
!-----------------------------------------------------------------------
      if (param(1) /= 4) call ERREUR (146)
      select case (nstoc(1))
      case (1)                                        ! CHERCHEC
         ns=nstoc(2)
         if (param(2) > 1) call ERREUR (147)
         call CHERCHEC(ns,tcon)
         if (param(3) == 3) then
            cte(ncte)=tcon
         else
            if (param(3) == 5) then
               call AJOUTC(nresul(3),tcon)
            else
               call ERREUR(147)
            end if
         end if
      case (2)                                        ! CONSER
         if (param(2) /= 3) call erreur (147)
         c=cstoc(2)
         nom=nomcte(nstoc(2))
         call KILLC(nom)
         call AJOUTN(nom)
         call QUELNUMERO(nom,n)
         call CONSER(c,n)
      case (3:4)
         if (nstoc(1)==3) xcos=.TRUE.                 ! HANSENC
         if (nstoc(1)==4) xcos=.FALSE.                ! HANSENS
         if (param(2) /= 1.AND.param(2) /= 5 ) call erreur(147)
         if (param(2) == 1) ns=nstoc(2)
         if (param(2) == 5) then
            call AJOUTN(nresul(2))
         ns=numser(nser)
         end if
         if (param(3) /= 2 .AND. param(3) /= 3) call erreur (147)
         nw=cstoc(3)
         if (param(4) /= 2 .AND. param(4) /= 3) call erreur (147)
         mw=cstoc(4)
         iordre=maxop
         call DEUXCORPS(nw,mw,xcos,iordre,ns)
      case (5:6)
         if (nstoc(1)==5) xcos=.TRUE.                 ! SINGLEC
         if (nstoc(1)==6) xcos=.FALSE.                ! SINGLES
         if (param(2) /= 1.AND.param(2) /= 5 ) call erreur(147)
         if (param(2) == 1) ns=nstoc(2)
         if (param(2) == 5) then
            call AJOUTN(nresul(2))
         ns=numser(nser)
         end if
         if (param(3) /= 2 .AND. param(3) /= 3) call erreur (147)
         nw=cstoc(3)
         if (param(4) /= 2 .AND. param(4) /= 3) call erreur (147)
         mw=cstoc(4)
         iordre=maxop
         call SINGLE(nw,mw,xcos,ns)
      case (7:8)
         if (param(2) > 1) call ERREUR (147)
         n1=nstoc(2)
         if (param(3) /= 1.AND.param(3) /= 5 ) call erreur(147)
         if (param(3) == 1) n2=nstoc(3)
         if (param(3) == 5) then
            call AJOUTN(nresul(3))
         n2=numser(nser)
         end if
         if (nstoc(1)==7) call DERARG(n1,n2)           ! DERARG
         if (nstoc(1)==8) call DERLIT(n1,n2)           ! DERLIT
      case (9)
         if (param(2) > 1) call ERREUR (147)
         n=nstoc(2)
         if (param(3) /= 1.AND.param(3) /= 5 ) call ERREUR(147)
         if (param(3) == 1) ns=nstoc(3)
         if (param(3) == 5) then
            call AJOUTN(nresul(3))
         ns=numser(nser)
         end if
         if (param(4) /= 1.AND.param(4) /= 5 ) call ERREUR(147)
         if (param(4) == 1) nc=nstoc(4)
         if (param(4) == 5) then
            call AJOUTN(nresul(4))
         nc=numser(nser)
         end if
         if (param(5) /= 2 .AND. param(5) /= 3) call ERREUR (147)
         norder=cstoc(5)
         call SINCOS(n,ns,nc,norder)                  ! SINCOS
      case default
         print *,'Appel de routine utilisateur:'
         print *,'Parametres: ',PARAM
         print *,'Numeros: ',NSTOC
         print *,'Constantes: ',CSTOC
         call ERREUR(148)
      end select
      END subroutine EXPLORE

!---------------------------------------------------------------------

!                             S I C O N S T

!---------------------------------------------------------------------
!     Test d'une chaine derriere le signe (=). Si cette chaine ne contient
!     que des constantes XTOT=.TRUE., sinon  XTOT=.FALSE.
!---------------------------------------------------------------------
      subroutine SICONST(chaine)
      use PARTLIT
      implicit real*8      (a-h,o-z)

      character*1          :: car
      character*20         :: nom
      character*100        :: chainbis
      character*(*)        :: chaine
!-----------------------------------------------------------------------
      chainbis=chaine
      xtot=.TRUE.
      nx=LEN_TRIM(chainbis)
      if (nx == 0) call ERREUR(140)
      do while (nx > 0)
         car=chainbis(1:1)
         select case (car)
         case ('$')
            xtot=.FALSE.
            RETURN
         case ('A':'Z')
            call SYMBOLE(chainbis,nom,nx)
            call TRIC(nom)
            if (xcte) goto 10
            call TRIN(nom)
            if (xser) then
               xtot=.FALSE.
               RETURN
            else
               call ERREUR(112)
            end if
         case ('0':'9','.')
            call NOMBRE(chainbis,c,nx)
         case ('+','-','*','/')
            nx=1
         case default
            call ERREUR(114)
         end select
!     Recherche d'un separateur (,) dans le cas d'une liste
10       continue
         chainbis=chainbis(nx+1:)
         nx=LEN_TRIM(chainbis)
      end do
      END subroutine SICONST


!***********************************************************************
!***********************************************************************

!     ---- G E S T I O N   D E S   T A B L E A U X  D E   N O M S ----

!***********************************************************************
!***********************************************************************

!----------------------------------------------------------------------

! [G]                         I N I T A B

!----------------------------------------------------------------------
!     Remise a zero de NOMSER(*) et NOMCTE(*) (caractere maximum) et
!     chargement de NDISMEM(*), NDISTEMP(*), des pointeurs de pile et
!     NMEM(*). Pour la signification des noms de tableaux cf. PUSH.
!----------------------------------------------------------------------
      subroutine INITAB
      use PARTLIT
      implicit real*8 (a-h,o-z)

!----------------------------------------------------------------------
      ntser=0
      ntcte=0
!     RAZ des noms/numeros de series
      do i=1,MXSER
      nomser(i)=BIGNOM
      numser(i)=0
      end do
!     RAZ des numeros de series en memoire
      do i=1,NBT
      ndismem(i)=NBT-i+1
      nmem(i)=0
      end do
      ntopmem=NBT
      ntmem=0
!     RAZ des numeros de series sur le tampon disque
      do i=1,MXTAMP
      ndistamp(i)=MXSER-i+1
      end do
      ntoptamp=MXTAMP
      do i=1,MXCTE
      nomcte(i)=BIGNOM
      cte(i)=0
      end do
!!!!! [G]      FREQ=0 ; PHASE=0
      END subroutine INITAB

!----------------------------------------------------------------------

! [G]                         T R I N

!----------------------------------------------------------------------
!     Tri dans la table NOMSER(*) par ordre alphanumerique croissant.
!     Entree: NOM a rechercher dans NOMSER(*);
!     Sortie: NSER et XSER. Si le NOM existe XSER=.TRUE., NSER est
!     le rang de NOM dans NOMSER(*). Dans le cas contraire XSER=.FALSE.,
!     NSER est l'adresse de NOM dans NOMSER pour un interclassement.
!----------------------------------------------------------------------
      subroutine TRIN(nom)
      use PARTLIT
      implicit real*8 (a-h,o-z)

      character*20                        :: nom,zw

!----------------------------------------------------------------------
      nser=milser
      xser=.TRUE.
      do i=idebser,ISTOP
      zw=nomser(nser)
      if (nom == zw) RETURN
      if (nom > zw) then
         nser=nser+NPUIS2(i)
         l=1
      else
         nser=nser-NPUIS2(i)
         l=0
      endif
      end do
      nser=nser+l
      xser=.FALSE.
      END subroutine TRIN

!----------------------------------------------------------------------

! [G]                         K I L L N

!----------------------------------------------------------------------
!     Le nom donne en entree (NOM) est elimine de la table NOMSER(*)
!     s'il existe, et la table est comprimee (NTSER=NTSER-1);
!     sinon la table reste inchangee.
!----------------------------------------------------------------------
      subroutine KILLN(nom)
      use PARTLIT
      implicit real*8 (a-h,o-z)


      character*20                        :: nom

!----------------------------------------------------------------------
      if (ntser == 0) RETURN
      call TRIN(nom)
      if (xser) then
      call PUSH(numser(nser))
      if (nser <  ntser) then
         do n=nser,ntser-1
            nomser(n)=nomser(n+1)
            numser(n)=numser(n+1)
         end do
      end if
      nomser(ntser)=BIGNOM
      numser(ntser)=0
      ntser=ntser-1
      else
         print *, 'KILLN: La serie de nom: ',nom,' n''existe pas'
      RETURN
      endif
      END subroutine KILLN

!----------------------------------------------------------------------

! [G]                         A J O U T N

!----------------------------------------------------------------------
!     Le nom donne en entree (NOM) est intercale dans la table NOMSER(*).
!     S'il n'existe pas deja la table est etendue (NTSER=NTSER+1);
!     sinon la table reste inchangee.
!     A chaque nouvel ajout d'un tableau en memoire NMEM(*) est charge
!     du numero de serie.
!----------------------------------------------------------------------
      subroutine AJOUTN(nom)
      use PARTLIT
      implicit real*8 (a-h,o-z)

      character*20                        :: nom

!----------------------------------------------------------------------
      call TRIN(nom)
      if (xser) then
      print *, 'AJOUTN: La serie de nom: ',nom,' existe deja'
      else
      if (ntser >= MXSER) call ERREUR(150)
      if (nser <= ntser) then
         do n=ntser,nser,-1
            nomser(n+1)=nomser(n)
            numser(n+1)=numser(n)
         end do
      end if
      nomser(nser)=nom
      call POP(numser(nser))
      ntser=ntser+1
      endif
      END subroutine AJOUTN

!----------------------------------------------------------------------

! [G]                         P U S H

!----------------------------------------------------------------------
!     Subroutines de manipulation de la pile memoire+tampon.
!     NDISMEM: pile des numeros de series disponibles en memoire;
!     NDISTAMP: pile des numeros de series diponibles dans le tampon.
!     NMEM: Liste des numeros de tableaux en memoire au nombre de NTMEM.
!----------------------------------------------------------------------
      subroutine  PUSH(n)
      use PARTLIT
      implicit real*8 (a-h,o-z)

!----------------------------------------------------------------------
      if (n <= NBT) then
      ntopmem=ntopmem+1
      if (ntopmem > NBT) call ERREUR(151)
      ntmem=ntmem-1
      do i=1,ntmem
         if (nmem(i) == n) then
            do j=i,ntmem
            nmem(j)=nmem(j+1)
            end do
            exit
         end if
      end do
      ndismem(ntopmem)=n
      else
      ntoptamp=ntoptamp+1
      if (ntoptamp > MXTAMP) call ERREUR(151)
      ndistamp(ntoptamp)=n
      end if
      END subroutine PUSH

!----------------------------------------------------------------------

! [G]                         P O P

!----------------------------------------------------------------------
!     Pour la signification des variables cf. la subroutine PUSH.
!----------------------------------------------------------------------
      subroutine POP(n)
      use PARTLIT
      implicit real*8 (a-h,o-z)

!----------------------------------------------------------------------
      if (ntopmem > 0) then
      n=ndismem(ntopmem)
      ntmem=ntmem+1
      nmem(ntmem)=n
      ntopmem=ntopmem-1
      else
      if (ntoptamp <= 0) call ERREUR(152)
      n=ndistamp(ntoptamp)
      ntoptamp=ntoptamp-1
      end if
      END subroutine POP

!----------------------------------------------------------------------

! [G]                         T R I C

!----------------------------------------------------------------------
!     Tri dans la table NOMCTE(*) par ordre alphanumerique croissant.
!     Entree: NOM a rechercher dans NOMCTE(*);
!     Sortie: NCTE et XCTE. Si NOM existe XCTE=.TRUE., NCTE est le rang
!     de NOM dans NOMCTE(*). Dans le cas contraire XCTE=.FALSE., NCTE est
!     l'adresse de NOM dans NOMCTE(*) pour un interclassement.
!----------------------------------------------------------------------
      subroutine TRIC(nom)
      use PARTLIT
      implicit real*8 (a-h,o-z)

      character*20                        ::  nom,zw

!----------------------------------------------------------------------
      ncte=milcte
      xcte=.TRUE.
      do i=idebcte,ISTOP
      zw=nomcte(ncte)
      if (nom == zw) RETURN
      if (nom > zw) then
         ncte=ncte+NPUIS2(i)
         l=1
      else
         ncte=ncte-NPUIS2(i)
         l=0
      endif
      end do
      ncte=ncte+l
      xcte=.FALSE.
      END subroutine TRIC

!----------------------------------------------------------------------

! [G]                         K I L L C

!----------------------------------------------------------------------
!     Le nom donne en entree (NOM) est elimine de la table NOMCTE(*)
!     s'il existe, et la table est comprimee (NTCTE=NTCTE-1);
!     sinon la table reste inchangee.
!----------------------------------------------------------------------
      subroutine KILLC(nom)
      use PARTLIT
      implicit real*8 (a-h,o-z)

      character*20                        :: nom

!----------------------------------------------------------------------
      if (ntcte == 0) RETURN
      call TRIC(nom)
      if (xcte) then
      if (ncte <  ntcte) then
         do n=ncte,ntcte-1
            nomcte(n)=nomcte(n+1)
            cte(n)=cte(n+1)
         end do
      end if
      nomcte(ntcte)=BIGNOM
      cte(ntcte)=0
      ntcte=ntcte-1
      else
      print *, 'KILLC: La cte de nom: ',nom,' n''existe pas'
      RETURN
      endif
      END subroutine KILLC

!----------------------------------------------------------------------

! [G]                         A J O U T C

!----------------------------------------------------------------------
!     Le nom donne en entree (NOM) est intercale dans la table NOMCTE(*)
!     S'il n'existe pas deja la table est etendue (NTCTE=NTCTE+1);
!     sinon la table reste inchangee.
!----------------------------------------------------------------------
      subroutine AJOUTC(nom,ct)
      use PARTLIT
      implicit real*8 (a-h,o-z)

      character*20                        :: nom

!----------------------------------------------------------------------
      call TRIC(nom)
      if (xcte) then
      print *, 'AJOUTC: La cte de nom: ',nom,' existe deja'
      RETURN
      else
      if (ntcte >= MXCTE) call ERREUR(153)
      if (ncte <= ntcte) then
         do n=ntcte,ncte,-1
            nomcte(n+1)=nomcte(n)
            cte(n+1)=cte(n)
         end do
      end if
      nomcte(ncte)=nom
      cte(ncte)=ct
      ntcte=ntcte+1
      endif
      END subroutine AJOUTC

!----------------------------------------------------------------------

! [G]                         Q U E L N O M

!----------------------------------------------------------------------
!     Recherche du nom d'une serie dans le tableau NOMSER(*) a partir
!     de son numero N dans NUMSER(*). Entree: N; Sortie: NOM.
!----------------------------------------------------------------------
      subroutine QUELNOM(n,nom)
      use PARTLIT
      implicit real*8 (a-h,o-z)

      character*20                        :: nom

!----------------------------------------------------------------------
      if (n == 0) then
         nom='$'
         RETURN
      end if
      do m=1,ntser
      if (numser(m) == n) then
         nom=nomser(m)
         RETURN
      end if
      end do
      call ERREUR(154)
      END subroutine QUELNOM

!-----------------------------------------------------------------------

!                             Q U E L N U M E R O

!-----------------------------------------------------------------------
!      Recherche du numero d'une serie dans le tableau NUMSER() a
!      partir de son nom dans NOMSER(). Entree: NOM; Sortie: N.
!-----------------------------------------------------------------------
      subroutine QUELNUMERO(nom,n)
      use PARTLIT
      implicit real*8      (a-h,o-z)

      character*20                        :: nom
!-----------------------------------------------------------------------
      do m=1,ntser
         if (nomser(m) == nom) then
            n=numser(m)
            RETURN
         end if
      end do
      call ERREUR(155)
      END subroutine QUELNUMERO


!***********************************************************************
!***********************************************************************

!     ---------------- G E S T I O N  D I S Q U E ----------------

!***********************************************************************
!***********************************************************************

!----------------------------------------------------------------------

! [G]                         C A S H

!----------------------------------------------------------------------
!     La subroutine teste le rang du numero NIN de serie.
!     Si NIN est un numero de serie en memoire, la serie est inchangee.
!     Si NIN est un numero de serie sur le tampon disque, la serie
!     en memoire la plus ancienne est stockee sur le tampon.
!     Si XLEC=.TRUE., serie NIN est lue et remplace cette derniere en
!     memoire. Si XLEC=.FALSE., NIN n'est pas lue.
!     NOUT est le numero de la serie qui remplace NIN.
!----------------------------------------------------------------------
      subroutine CASH(nin,xlec,nout,interdit)
      use PARTLIT
      implicit real*8 (a-h,o-z)

      logical*1                           :: xlec

      character*20                        :: noma,nomb

!----------------------------------------------------------------------
      NLOG=NLOG3
      if (nin == 0) then
      nout=0
      RETURN
      endif
      if (nin <= NBT) then
      mm=ntmem-1
      do i=1,mm
         if (nmem(i) == nin) then
            do j=i,mm
            nmem(j)=nmem(j+1)
            end do
            exit
         end if
      end do
      nmem(ntmem)=nin
      nout=nin
      RETURN
      end if
      call QUELNOM(nin,noma)
      call TRIN(noma)
      nsa=nser
      if (ntmem <  NBT) then
      call POP(numser(nsa))
      call PUSH(nin)
      nb=numser(nsa)
      if (xlec) call LDISK(NLOG,nb,nin-NBT+1,noma)
      nout=nb
      else
      do i=1,ntmem
            nb=nmem(i)
            if (nb /= interdit) exit
         end do
      call QUELNOM(nb,nomb)
      call TRIN(nomb)
      nsb=nser
      numser(nsa)=nb
      if (xlec) then
         call POP(numser(nsb))
         nc=numser(nsb)
         call EDISK(NLOG,nb,nc-NBT+1,nomb)
         call LDISK(NLOG,nb,nin-NBT+1,noma)
         call PUSH(nin)
      else
         numser(nsb)=nin
         call EDISK(NLOG,nb,nin-NBT+1,nomb)
      end if
         nout=nb
         mm=ntmem-1
         do i=1,mm
            if (nmem(i) == nb) then
               do j=i,mm
                  nmem(j)=nmem(j+1)
               end do
               exit
            end if
         end do
         nmem(ntmem)=nb
      end if
      END subroutine CASH

!----------------------------------------------------------------------

!                          E D I S K

!----------------------------------------------------------------------
!     Ecriture d'une serie de numero N, a l'enregistrement NREC.
!----------------------------------------------------------------------
      subroutine EDISK(NLOG,n,ndisk,nom)
      use PARTLIT
      implicit real*8 (a-h,o-z)
      character*20                        :: nom

!----------------------------------------------------------------------
      nrec=3*(ndisk-1)+1
      write(NLOG3,rec=nrec) (t(i,n),i=1,MXDIMT)
      nrec=nrec+1
      write(NLOG3,rec=nrec) (jt(i,n),i=1,MXDIMT),(lt(i,n),i=1,MXDIMT)
      nrec=nrec+1
      write(NLOG3,rec=nrec) nom,ntt(n),trig(n)
      RETURN
      END subroutine EDISK

!----------------------------------------------------------------------

!                          L D I S K

!----------------------------------------------------------------------
!     Lecture d'une serie de numero N, a l'enregistrement NREC.
!----------------------------------------------------------------------
      subroutine LDISK(NLOG,n,ndisk,nom)
      use PARTLIT
      implicit real*8 (a-h,o-z)
      character*20                        :: nom,nomdisk

!----------------------------------------------------------------------
      nrec=3*(ndisk-1)+1
      read(NLOG3,rec=nrec) (t(i,n),i=1,MXDIMT)
      nrec=nrec+1
      read(NLOG3,rec=nrec) (jt(i,n),i=1,MXDIMT),(lt(i,n),i=1,MXDIMT)
      nrec=nrec+1
      read(NLOG3,rec=nrec) nomdisk,ntt(n),trig(n)
      if (nomdisk /= nom) then
         print *,'NOMDISK: ',nomdisk,' NOM: ',nom
         call ERREUR(160)
      end if
      RETURN
      END subroutine LDISK

!***********************************************************************
!***********************************************************************

!     ------------ R O U T I N E S  S P E C I A L E S ------------

!***********************************************************************
!***********************************************************************

!-----------------------------------------------------------------------

!                             C H E R C H E C

!-----------------------------------------------------------------------
      subroutine CHERCHEC(ns,tcon)
!-----------------------------------------------------------------------
!     Recherche du terme constant d'une serie dans un tableau T ou $
!     Entree : NS : numero du tableau
!     Sortie : TCON : valeur du terme constant
!-----------------------------------------------------------------------
      use PARTLIT
      implicit real*8      (a-h,o-z)

      logical                             :: xcomp
!-----------------------------------------------------------------------
      call CASH(ns,.TRUE.,n1,0)

      tcon=0
      xcomp=.TRUE.
      if (n1 == 0) then
         if (.NOT.trigr) call ERREUR(245)
         do i=1,ntr
            xcomp=xcomp.AND.(jr(i)==0).AND.(lr(i)==0)
            if (xcomp) then
               tcon=r(i)
               EXIT
            end if
         end do
      else
         if (.NOT.trig(n1)) call ERREUR(246)
         do i=1,ntt(n1)
            xcomp=xcomp.AND.(jt(i,n1)==0).AND.(lt(i,n1)==0)
            if (xcomp) then
               tcon=t(i,n1)
               EXIT
            end if
         end do
      end if

100   ns=n1
      END subroutine CHERCHEC

!-----------------------------------------------------------------------

!                             C O N S E R

!-----------------------------------------------------------------------
      subroutine CONSER(c,n)
!-----------------------------------------------------------------------
!     Transformation d'une constante C en une serie C*cos(0) de numero N.
!-----------------------------------------------------------------------
      use PARTLIT
      implicit real*8      (a-h,o-z)

!-----------------------------------------------------------------------
      call CASH(n,.FALSE.,n1,0)
      call RAZ(n1)
      trig(n1)=.TRUE.
      lt(1,n1)=0 ; jt(1,n1)=0
      t(1,n1)=c
      ntt(n1)=1

      n=n1
      END subroutine CONSER


!***********************************************************************
!***********************************************************************

!     ----------- G E S T I O N   D E S   E R R E U R S ----------

!***********************************************************************
!***********************************************************************

!----------------------------------------------------------------------

!                          E R R E U R

!----------------------------------------------------------------------
!     Routine de gestion des erreurs
!----------------------------------------------------------------------
      subroutine ERREUR(n)
      use PARTLIT
      implicit real*8 (a-h,o-z)

!---------------------------------------------------------------------
      write (*,'(a,i4,a)') ' --------- E R R E U R  No: ',n,'----------'

      select case (n)
      case (100)
      print *,'ZONE DE LECTURE: ',zrecbis
      case (105:149)
      print *,'INSTRUCTION: ',instrucbis
      end select

      select case (n)
      case (100)
      STOP 'Limite de commentaire incomplete'
!      case (105)
!     STOP 'Numero de parametre formel trop grand'
      case (106)
      STOP 'Nom de constante reservee inconnu'
      case (109)
      STOP  'Valeur de SIFRAC ou SIARGU incorrecte'
      case (110)
      STOP 'Cette instruction n''est pas une affectation'
      case (111)
      STOP 'La serie n''existe pas'
      case (112)
      STOP 'Le symbole cte/serie n''existe pas'
      case (113)
      STOP 'La commande n''existe pas'
      case (114)
      STOP 'Lecture d''un caractere inconnu'
!     case (115)
!     STOP 'Le numero de variable/parametre n''existe pas'
      case (116)
      STOP 'La variable/parametre est hors des limites'
!     case (117)
!     STOP 'Parametres des 2-corps absents'
!     case (118)
!     STOP 'Parametres de la saisie d''une serie simple absents'
      case (119)
      STOP 'Operation de substitution du temps impropre'
      case (120)
      STOP 'Operation impossible avec $'
      case (121)
      STOP '2 constantes dans un seul cycle'
      case (122)
      STOP 'La variable du membre droit n''est pas une serie'
      case (123)
      STOP 'Plus de 3 elements dans un cycle'
      case (124)
      STOP 'Operateur inconnu'
      case (125)
      STOP 'Cycle operatoire trop long'
      case (126)
      STOP 'Lecture d''un caractere inconnu'
!     case (127)
!     STOP 'Ligne trigonometrique absente'
      case (128)
      STOP 'La serie donnee n''existe pas sur disque'
      case (129)
      STOP 'Nom de fichier absent'
      case (130)
      STOP 'Affectation serie=cte impossible'
      case (131)
      STOP 'Affectation cte=serie impossible'
      case (132)
      STOP 'Operation impossible avec $'
      case (133)
      STOP 'erreur de numero de cycle'
      case (136)
      STOP 'Pas de debut de boucle dans une iteration'
      case (137)
      STOP 'Pas de fin de boucle | Iterations imbriquees'
      case (138)
      STOP 'Operation sur ctes: la constante n''existe pas'
      case (139)
      STOP 'Division de constante nulle'
      case (140)
      STOP 'Chaine vide'
      case (141)
      STOP 'Longueur de nom trop longue'
      case (142)
      STOP 'La chaine n''est pas un nombre'
      case (143)
      STOP 'Longueur de nombre trop longue'
      case (145)
      STOP 'Variable inconnue dans la liste (routine utilisateur)'
      case (146)
      STOP 'Nom de routine utilisateur absent en 1er parametre'
      case (147)
      STOP 'Appel de la routine utilisateur incorrect'
      case (148)
      STOP 'Subroutine utilisateur inexistante dans EXPLORE'
      case (149)
      STOP 'Subroutine utilisateur inexistante dans FILSUB'
      case (150)
      STOP 'La table des noms est saturee'
      case (151)
      STOP 'Sur-depassement de pile (PUSH)'
      case (152)
      STOP 'Sous-depassement de pile (POP)'
      case (153)
      STOP 'La table des constantes est saturee'
      case (154)
      STOP 'Il n''exite pas de serie avec ce numero'
      case (155)
      STOP 'Il n''existe pas de numero avec ce nom de serie'
      case (160)
      STOP 'Il n''existe pas de serie sur disque avec ce nom'
      case (170)
      STOP  'Operation sur les constantes illicite'
      case (171)
      STOP 'Il n''existe pas de symbole avec ce nom dans l''impression'
      case (200)
      STOP 'Nombre d''arguments trigonometriques trop grand'
      case (201)
      STOP 'Nombre de parametres litteraux trop grand'
      case (202)
      STOP 'Bornes litterales trop grandes'
      case (203)
      STOP 'Ordre du n-uplet trop grand [TASSE]'
      case (204)
      STOP 'Ordre du n-uplet trop grand [DETASSE]'
      case (205)
      STOP 'Adressage incorrect [TRIX]'
      case (206)
      STOP 'Accumulateur $ sature'
      case (207)
      STOP 'Adressage incorrect [DECAL]'
      case (208)
      STOP 'Adressage incorrect [TRIFI]'
      case (209)
      STOP 'Table T(*,*) saturee'
      case (210)
      STOP 'Adressage incorrect [DECLIT]'
      case (211)
      print  *, 'Ordre de developpement trop grand dans (SINCOS)'
      case (221)
      STOP 'PURGE de $ non autorisee'
      case (222)
      STOP 'Lignes trigonometriques incompatibles [MULSER]'
      case (223)
      STOP 'Lignes trigonometriques incompatibles [SOMSER]'
      case (224)
      STOP 'Somme d''une constante a une serie en sinus interdite'
      case (225)
      STOP 'Boucle infinie [FRACONT]'
      case (226)
      STOP 'Coefficient de Hansen d''ordre trop grand'
      case (227)
      STOP 'erreur dans le calcul de la fonction hypergeometrique'
      case (230)
      STOP 'INTEGRATION: Serie de Fourier de terme constant /=0'
      case (240)
      STOP 'Derivation de $ interdite dans DERIVE'
      case (241)
      STOP 'Integration de $ interdite dans PRIMIT'
      case (242)
      STOP 'Derivation partielle de $ interdite dans DERARG'
      case (243)
      STOP 'Derivation partielle de $ interdite dans DERLIT'
      case (244)
      STOP 'Substitution du temps dans $ interdite dans SUBTIME'
      case (245)
      STOP 'Il n''y a pas de terme constant dans une serie en sinus'
      case (246)
      STOP 'Il n''y a pas de terme constant dans une serie en sinus'
      case (247)
      STOP 'Tranformation de $ en une serie GREGOIRE interdite'
      case (250)
      STOP 'Troncature operationnelle (MAXOP) > Ordre max. (MAXORD)'
      case (251)
      STOP 'Liste erronee dans la recherche de coefficients/arguments'
      case default
      STOP 'Numero d''erreur inconnu'
      end select
      END subroutine ERREUR

!-----------------------------------------------------------------------

!                             S I N C O S

!-----------------------------------------------------------------------
      subroutine SINCOS(n,ns,nc,norder)
!-----------------------------------------------------------------------
!     Calcul du SINUS ou du COSINUS d'une serie. On suppose que la serie
!     a la forme suivante S = s_0 + Ds. s_0 est la constante de la serie,
!     Ds son developpement. Ds est une petite quantite.
!     s_0 et Ds sont exprimes en radian. On forme
!     sin(s_0+Ds) = sin(s_0)cos(Ds)+cos(s_0)sin(Ds), et,
!     cos(s_0+Ds) = cos(s_0)cos(Ds)-sin(s_0)sin(Ds).
!     On developpe en serie de Taylor jusqu'a l'ordre NORDER, puissance
!     maximum de Ds, soit
!     sin(Ds) = Ds - 1/3!Ds**3 + 1/5!Ds**5 - ...
!     cos(Ds) = 1  - 1/2!Ds**2 + 1/4!D**4 +
!     On impose NORDER <= NORDERMAX
!     La constante s_0 est extraite de la serie initiale Ds
!
!     n : numero de la serie DS
!     ns, nc : numeros des series SINUS et COSINUS respectivement.
!-----------------------------------------------------------------------
      use PARTLIT
      implicit real*8      (a-h,o-z)
!     Puissance max. dans les developpement de Taylor SINUS/COSINUS de series
      integer, parameter   :: NORDERMAX=7
      logical*1            :: xsinus,xjump
      character*20         :: nom1='xxxxxxxxxxxxxxxxxxxx'
      character*20         :: nom2='yyyyyyyyyyyyyyyyyyyy'
!-----------------------------------------------------------------------

      if (NORDER > NORDERMAX) call ERREUR (211)
      call AJOUTN(nom1)
      n1=numser(nser)
      call AJOUTN(nom2)
      n2=numser(nser)
      if (trig(n)) then
         call CHERCHEC(n,s_0)
         cos_0=cos(s_0) ; sin_0=sin(s_0)
         call SOMCTE(n,ns,0,-s_0)
         xjump=.TRUE.
      else
         call TRANS(n,ns)
         xjump=.FALSE.
      end if

      UN=1.D0
      call CONSER(UN,nc)
      factr=UN
      isign=-1
      xsinus=.FALSE.
      call TRANS(ns,n2)
      call TRANS(ns,n1)
      do it=2,norder
         factr=UN/it
         f=isign*factr
         call MULSER(n2,n1,n1,0,f)
         if (xsinus) then
            call SOMSER(ns,ns,1,UN)
            isign=-isign
         else
            call SOMSER(nc,nc,1,UN)
            if (it == 2) then
               isign=1
            else
               isign=-isign
            end if
         end if
         xsinus=.NOT.xsinus
      end do
      if (xjump) then
         call SOMSER(nc,0,0,cos_0)
         call SOMSER(ns,n1,1,-sin_0)
         call SOMSER(nc,0,0,sin_0)
         call SOMSER(ns,ns,1,cos_0)
         call TRANS(n1,nc)
      end if

      call KILLN(nom1)
      call KILLN(nom2)

      END subroutine SINCOS
